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Abstract 

Previous studies have shown that the Hidden Eocal Symmetry (HES) Model, sup¬ 
plied with appropriate symmetry breaking mechanisms, provides an Effective Eagrangian 
(BHES) which encompasses a large number of processes within a unified framework. 
Based on it, a global fit procedure allows for a simultaneous description of the e"''e“ an¬ 
nihilation into 6 final states - 7r'''7r“, 7r°7, 77, 7r+7r“7r*^, K^K~, KlKs - and includes 
the dipion spectrum in the r decay and some more light meson decay partial widths. The 
contribution to the muon anomalous magnetic moment of these annihilation channels 
over the range of validity of the HES model (up to 1.05 GeV) is found much improved in 
comparison to the standard approach of integrating the measured spectra directly. How¬ 
ever, because most spectra for the annihilation process e''"e“ —)• 7r'''7r“ undergo overall 
scale uncertainties which dominate the other sources, one may suspect some bias in the 
dipion contribution to which could question the reliability of the global fit method. 
However, an iterated global fit algorithm, shown to lead to unbiased results by a Monte 
Carlo study, is defined and applied succesfully to the e'''e“ tt+tt” data samples from 
CMD2, SND, KEOE, BaBar and BESSIII. The iterated fit solution is shown to further 
improve the prediction for a^, which we find to deviate from its experimental value above 
the 4 c 7 level. The contribution to of the 7r'''7r“ intermediate state up to 1.05 GeV has an 
uncertainty about 3 times smaller than the corresponding usual estimate. Therefore, global 
fit techniques are shown to work and lead to improved unbiased results. 
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1 Introduction 


As is well known, the Standard Model is the gauge theory whieh eovers the realm of weak, 
eleetromagnetie and strong interaetions among quarks, leptons and the various gauge bosons 
(gluons, photons, W^, Z^). In energy regions where perturbative methods apply, the Standard 
Model (SM) allows to yield preeise estimates for several physical effects, sometimes with ac¬ 
curacies of the order of a few 10“^^. In contrast, in energy regions where the non-perturbative 
regime of QCD is involved, getting similar precision may become challenging. This is the case 
for the low energy part of the photon hadronic vacuum polarization (HVP); this HVP plays a 
crucial role in determining the theoretical value for the muon anomalous moment a^, one of 
the best measured particle properties. 

Fortunately, getting precise estimates in the low energy hadron SM sector is not completely 
out of reach as exemplified by the Chiral Perturbation Theory (ChPT)[|Tll3 which is rigorously 
the low energy limit of QCD, valid up to 400 500 MeV but lets the resonance region outside 

its scope. Lattice QCD (LQCD) is also a promising method under rapid development which 
already allows to perform precise computations at low (and very low) energies Q . Interesting 
LQCD estimates for the HVP’s of the three leptons have already been produced [jH [51 which 
clearly show that LQCD reaches results in accord with expectations; this is especially striking 
for with, however, still unsatisfactory uncertainties Q. 

So, much progress remains to be done before LQCD evaluations can compete with the 
accuracy of the experimental measurements already available [|S 13 or, a fortiori, with those 
expected in a near future at Fermilab [[8l|9l or, slightly later, at J-PARC ifTOl . Since lattice 
QCD is intrinsically an Euclidean approach, it is intrinsically unable to account for the existing 
rich amount of low energy hadronic data in the non-perturbative time-like region, i.e. from 
thresholds to 2 3 GeV. Therefore, other methods, able to encompass large fractions of the 

physics from this important energy region, are valuable. 

A natural approach to this issue is provided by Effective Lagrangians which cover the reso¬ 
nance region. Such Effective Eagrangians should be constructed so as to preserve the symmetry 
properties of QCD as already done by standard ChPT, however only valid up to the r] mass re¬ 
gion. As it includes meson resonances, the Resonance Chiral Perturbation Theory (R^PT) ifTTI 
is an appropriate framework to study annihilations from their respective thresholds up to 
the intermediate energy region. 

It has been proven lfT3 that the coupling constants occuring at order in ChPT are sat¬ 
urated by low lying meson resonances of various kinds as soon as they can contribute. This 
emphasizes the role of the fundamental vector meson nonet (V) and confirms the relevance of 
the Vector Meson Dominance (VMD) concept in low energy physics. 

On the other hand, it has been proven [fTTI that the Hidden Eocal Symmetry (HES) Model 
[[T3]l and R^PT are equivalent provided consistency with the QCD asymptotic behavior is in¬ 
corporated. It thus follows that the HES model is also a motivated and constraining QCD 
rooted framework. As the original HES Model only deals with the lowest mass resonances, it 
provides a framework for the annihilations naturally bounded by the f mass region - i.e. 
up to ~ 1.05 GeV. 

The non-anomalous [[T4ll and anomalous ifTSil sectors of the HES Model open a wide scope 
and can deal with a large corpus of physics processes in a unified way. However, as such, 
HES cannot precisely reach the numerical precision requested by the wide ensemble of high 
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statistics data samples collected by several sophisticated experiments on several annihilation 
channels. In order to achieve such a program, the HLS Lagrangian must be supplied with 
appropriate symmetry breaking mechanisms not present in its original formulation [[T3l . 

This was soon recognized by the HLS Model authors who first proposed the mechanism 
to break SU(3) symmetry IfTbl named BKY according to its author names. Its success was 
illustrated by several phenomenological studies based on the BKY breaking scheme IfTTl ITSl 
[T^ . It was also soon extended to SU(2)/Isospin symmetry breaking llIOll . However, in order 
to account simultaneously for all the radiative decays of the light flavor mesons, the additional 
step of breaking the nonet symmetry for light pseudoscalar mesons was required; based on the 
heuristic formulation of the VP'y couplings by O’Donnell [l2n which includes nonet symmetry 
breaking in the pseudoscalar (P) sector in a specific way, a global and successful! account of 
all VP'y and P 77 couplings has been reached ll^ . The BKY SU(3) breaking and this nonet 
symmetry breaking included within the HLS Model was shown [|2^ to meet the requirements of 
Extended Chiral Perturbation Theory [|2^ 1^ . Finally, introducing the physical vector meson 
fields as the eigenstates of the loop modified vector meson mass matrix provided a mixing 
scheme of the — cc — 0 system which together with the E — 7 loop transitions implied by 
the HLS model at one loop0 leads to a satisfactory solution [l27ll of the long-standing r — e+e“ 
puzzle [|28l|2^ [3011311. 

Therefore, the approach just sketched is a global framework aiming at accounting for the 
largest possible ensemble of data spectra collected in the largest possible number of low energy 
physics channels. As this global model is an Effective Eagrangian constructed from the (P 
and V) fields relevant in the low energy regime of QCD and because it is consistent with the 
symmetries of QCD, one naturally expects their low energy results to be consistent with the 
Standard Model. 

It was then shown, that the Effective Eagrangian constructed from the original HES Model 
supplemented with the breaking schemes listed above was able to provide a satisfactory simul¬ 
taneous description of the e+e“ annihilations into the 7 r’'' 7 r“, 7 r° 7 , 77 , 7 r+ 7 r“ 7 r° final states and 
of the dipion spectrum in the decay of the r lepton [|^ |33l. This tended to indicate that the 
T — e+e“ puzzle just referred was related to an incomplete incorporation of isospin symmetry 
breaking effects within models. 

Slightly extending these breaking schemes, one is led to the Broken HES (BHES) Model 
[|34l which provides a fully consistent picture of all examined annihilation cross section^ 
the T dipion spectrum and, additionally, some light meson decay information with a limited 
number of free parameters to be extracted from data. An interesting outcome of the BHES 
based fit framework was a novel evaluation of the dominating low energy piece of the HVP, 
leading to an improved estimate of the muon anomalous magnetic moment at more than 4cr 
from its measured valued [jbllTl. 

Introducing the dipion spectra collected in the ISR mode confirmed that the muon g — 2 
departs from expectation by more than Aa [|35]| . One should note that the high statistics ISR 

'See also li^ where the role of the — 7 mixing is especially emphasized. 

^Specifically the 6 e+e” annihilation channels to tt+tt”, 7r°7, 777, 7r+7r“7r°, K'^K~, K'^K^, each from its 
threshold up to 1.05 GeV, i.e. including the (j) signal region. 

^One should note that the BHLS evaluation for the muon HVP is the closest to the central value preferred by 
the Lattice QCD study ID . 
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dipion spectra recently published by the KLOE [[^ [J7l [38]|, BaBar [[^ l40l and BESSIII 
ll4ll Collaborations are strongly dominated by overall scale {i.e. normalization) uncertainties; 
additionally the KEOE and BaBar normalization uncertainties are energy dependent. However, 
sizeable overall scale uncertainties raise an important issue related with their possibly biasing 
the physics quantity values extracted from their spectra. This issue has been identified in the 
reference work of G. D’Agostini where a very simple case is proposed which illustrates 
that biasing effects can be dramatic^. Of course, for a key quantity like the muon g — 2, the 
problem should be explored and possible biases identified and fixed. The way out is already 
mentioned in [|4^ and further emphasized in other studies [|4?1 l46l @41 ; the exact solution 
exhibits a delicate issue as the removal of the bias on some quantity supposes to know its exact 
value. Nevertheless, as already suggested in ll42l and emphasized in ll44l . iterative methods can 
be defined and are expected to be bias free; this has been applied successfully to the derivation 
of parton density functions in ll47l . 

The present work mostly aims at reexamining the results provided in [l34l [35l concerning 
the muon HVP using an appropriately defined iterative fit method adapted to the dealing with 
form factors or cross sections in such a way that fit results and derived quantities - like the 
HVP, but not only - could be ascertained to be bias free. In this way, one can positively answer 
the question raised in the title of this study at the methodological level. 

The real issue of the physics model dependence can only be answered by having at disposal 
results derived from several independent model frameworks, all successfully (undoubtfully) 
accounting for the largest possible corpus of data. Indeed, the physics correlations relating the 
different physics processes encompassed within a given framework cannot easily accomodate 
a model-independent approach. Moreover, several issues within the global fit approach are 
related with the formulation of isospin symmetry breaking (IB) which can hardly be made 
model independent, especially in a global framework. 

The paper is organized as follows. In Section |2l one briefly reminds the concern of using 
Effective Eagrangian global frameworks in order to strengthen the constraints on the parameters 
to be derived from global fits. As our HES Eagrangian framework has a range limited upward 
to 1.05 GeV, the brief Section |3] reminds how the full HVP is derived from fit results and from 
additional information. 

Section @]is, actually, the center piece of the present paper as its purpose is to define the fit 
method when one should deal with samples affected by strong overall scale uncertainties. This 
firstly turns out to precisely define the functions to be minimized, depending on the specific 
properties of the spectra considered and, secondly, to set up and justify the iterative procedure 
we propose^. Subsection 14. 2l puts special emphasis on the specific function associated with 
samples affected by overall scale uncertainties besides a more usual experimental error matrix. 
The iterative fit procedure to deal with biases is formulated therein. 

Most of the ISR data samples exhibit s-dependent overall scale uncertainties, which are 
certainly a novel feature in our field; Subsection l4.Si delines an appropriate function suitable 
for such a case. Einally, Subsection 14.41 reports on the main features of the iterative global fit 

"'The issue raised by G. D’Agostini in this paper has also been met formerly in the context of Nuclear Physics 
where it is referred to as the ’’Peelle’s Pertinent Puzzle” (PPP) ll 43 ll which is examined thoroughly in HTll . 

^After completion of this work, we found that Bsl applies a method similar to ours to derive unbiased parton 
density functions from various kinds of measured spectra. 
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method when fitting sets of data samples eontaining samples with overall seale uneertainties 
of various magnitudes eompared to statistieal errors. The eonelusions reported here rely on a 
Monte Carlo study outlined and illustrated in Appendix lAl 

Seetion [5] reminds the data samples used within the BHLS proeedure and reports for a 
(minor) eorreetion affecting the amplitudes for the annihilation channels 7 r °7 and 77 . Section 
reports on the updated results of the fits performed using only the scan data and discarding all 
ISR data samples; the effects of the iterative method is illustrated here and it is shown that the 
needed number of iterations in the global fit procedure does not exceed 1. The more general 
running is the subject of Section |7] where updated results are given to correct for coding bugs 
affecting some of the numbers given in our [f34l |35l . The properties of the recently published 
KLOE12 [|^ and BESSIII iHTl data samples are examined. The evaluation of the muon g — 2 
based on the iterated fits of various combinations of data samples is the subject of Section [ 8 l 
where the HVP slope at s = 0 is also computed within BHES and compared to its value directly 
derived from experimental data. Einally, Section|9]is devoted to conclusions and remarks. 


2 Effective Lagrangian Frameworks And Global Fits 

As reminded in the Introduction, it is a common approach to rely on the Effective Ea- 
grangian (EE) method to cover the low energy region where QCD exhibits its non-perturbative 
regime and where the quark and gluon degrees of freedom are replaced by hadron fields. Each 
EE of practical use generally depends on parameters originating from the starting Eagrangians 
(like the pion decay constant or the universal vector coupling g) and on parameters gen¬ 
erated by the unavoidable symmetry breaking effects (like quark mass differences); all such 
parameters are determined from data with various precisions. 

Needless to say that any (broken) Effective Eagrangian provides amplitudes expected to 
account simultaneously for several different processes. This has a trivial consequence which, 
nevertheless, deserves to be stressed : All the Effective Eagrangians predict physics correlations 
among the different physical processes they can encompass : EL = {ifj, 1 = 1, • • • p}. 

Therefore, having plugged from start the physics correlations inside the (broken) Eagrangian, 
the amplitudes derived herefrom should allow for a global, simultaneous and constrained fit of 
all available data samples covering all the channels in EL. Provided the global fit is clearly 
successful, the parameter central values and uncertainties returned can be considered as the 
optimal values accounting for all the processes in Ei simultaneously. Therefore, one can con¬ 
sider that the fit information - parameter central values and error covariance matrix - exhausts 
the experimental information contained in the data samples covering all the processes in EL. 

Prom now on, one specializes to the Broken HES (BHES) model as defined and used in 
dsn. All data samples used in the global fit procedure defined in this paper have already 
been listed and analyzed in this Referencel^; this will not be repeated here. As for the vr+vr” 
annihilation final state, which is a central piece of HVP studies, this Reference dealt with only 
the available scan data which are dominated by the samples from CMD2 [l52l [5^ and SND 

® Concerning the non-7r+7r“ channels, all existing data samples collected in scan mode at Novosibirsk are 
considered. The r data included in the global fit procedure are the samples collected by ALEPH 1491 . CLEO l50l 
and Belle ED. 
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If54l . The samples eollected in the ISR mode by Babar lf55l as well as the former KLOE data 
samples (KLOE08 [[^ and KEOEIO (371) have been considered in (351 . Preliminary results 
including also the most recent KEOE sample (KEOE12) (3^ have been given in (5^1571 . The 
BESSIII spectrum iHTl . published by mid of 2015, is also included within our analysis. 

3 Estimating the Muon Non-Perturbative HVP 

The issue raised in this paper is whether Effective Eagrangian methods really improve the 
evaluation of the dominating non-perturbative part of the HVP (341 [35l compared to a direct 
integration of experimental data (see (2^ [5^ [591 for instance). As we are working within the 
original HES framework (131 . what is discussed is the HVP fraction associated with the tt+tt”, 
7 r°7, 77 , 7 r+ 7 r“ 7 r°, K^K~, intermediate states - covered by BHES - up to ~ 1.05 GeV; 
this represents more than 80% of the total EO-HVP . 

Basically, the leading order (EO) non-perturbative QCD contribution to the muon HVP is 
estimated separately for each intermediate hadronic state Hi via : 



( 1 ) 


and the total non-pertubative HVP component is the sum of all the possible a^(ifi). The 
function K{s) in Eq. dU) is a known kernel (311 enhancing the threshold regions (s^i) for any 
channel Hi and <JHi{s) is the undressed cross section 0 for the e+e“ —)■ Hi annhilation; Scut 
is an energy limit above which perturbative expansions are supposed to become valid. BHES 
permits to evaluate the 6 integrals z = 1, • • ■ 6 } up to ~ 1.05 GeV. As the energy 

interval [s^,Scut\ contribution to a^{Hi) is beyond the BHES energy range of validity, it is 
estimated using customary methods (like those defined in [|58l[59l[^, for instance), as also 
the full contributions of the channels outside the present BHES scope, like the 4 pion final 
states. As already stated, these pieces represent altogether about 20% of the muon EO-HVP 
contribution to a^. 

As can be checked by looking at the cross section formulae given in (341 . most parameters 
to be fitted appear simultaneously in the 6 different cross sections {cx//. (s), z = 1 , • • ■ 6 } and 
each annihilation channel if* comes in with several experimental data sample^. Therefore, for 
instance, the data samples covering any of the vr^y, 77 , K^K~, annihilation 

channels play as additional constraints on the tt+tt" cross section and are treated on the same 
footing than the annihilation data themselves. On the other hand, the constraints carried 
by the dipion r decay spectrum data (49l [50l ISTI influence the fit and allow to reduce the 
BHES parameter uncertainties in a consistent wajo This explains why the global fit method is 
expected to improve each a^{Hi) contribution compared to more traditional methods - those 
from (2^ [5^ [59l for instance - as these ignore the inter-channel correlations revealed by the 
BHES Effective Eagrangian and validated by satisfactory global fits. Of course, inter-channel 

^Final state radiation (FSR) effects also contribute and are estimated as in M- 

® An experimental data sample is defined as the measured spectrum m and all the uncertainties which affect it. 

®So also do the decay partial widths of the form P —^ 77 and V P7 (or 77' — coj) extracted from the 
Review of Particle Properties (RPP) libTl and implemented within BHLS. 
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correlations are a general feature of Effeetive Lagrangians, and not partieular for the BHLS 
implementation. 

As any method, the BHLS based global fit method earries speeifie systematies whieh have 
been examined in great details in If35l . It is worth remarking, to avoid ambiguities, that the 
isospin breaking effeets speeifie of the r dipion speetra are introdueed in the dipion speetrum 
[1^ as eommonly done in the literature |63l |Ml |651 |Ml |671 |Ml IMl (see also f[2^1 : they 
are totally independent of the isospin breaking sehemes involved in the BHLS Lagrangian and, 
aetually, eome supplementing these [f35l . 

4 Can One Trust Global Fit Results ? 

The global fit method previously used in If34ll35ll defines a so-ealled VMD strategy whieh 
ean be phrased in the following way : 

• 1/ If the physies eorrelations predieted by a given Effeetive Lagrangian Model are sup¬ 
ported by the experimental data they eneompass, they ean be eonsidered as exaet at the 
accuracy level reported for the data. 

• 2/ Whenever the deseription - global fit - provided by a given Effeetive Lagrangian is 
satisfaetory, the model eross seetions, the fit parameter values and the parameter error 
eovarianee matrix exhaust reliably the physies information eontained in the fitted data 
samples. 

In the present ease where the BHLS model is eoneerned, and foeusing on the muon LO- 
HVP, Statement # 2 means that the improvements for the 6 aeeessibles a^{Hi) derived from Eq. 
([U) by integrating from to 1.05 GeV/e are legitimately valid and eoneeptually supported. 

On the other hand. Statement # 1 does not mean that the importanee of the word ’’Effeetive” 
is forgotten, as elear from the italie sentenee it earries : Its validity might have to be revised if 
the experimental eontext evolves towards a degraded aeeount of the datc0. 

Obviously, a VMD strategy heavily relies on the statistieal methods used to analyze and 
fit the data; thus, one should aseertain that all aspeets of the data handling are taken into ae- 
eount as they should. In partieular, all features of the experimental uneertainties should be 
implemented eanonieally within the minimized global and in the fitting proeedure. Indeed, 
as remarked in [|451 ITOl . ineorreet fit results are more frequently due to an ineorreet dealing 
with the experimental errors (and eorrelations) rather than to the minimization proeedure itself. 
Therefore, speeial eare is requested in dealing with experimental uneertainties and in ehoosing 
the appropriate expression adapted to eaeh data sample. 

It is the purpose of this Seetion to address this issue and eheek whether the proeedure 
defined in [l34l |351 fulfills this statement; this will lead us to eomplement the fitting proeedure 
by an iterative method. 

'^However, if an ensemble of data is internally conflicting within a given Effective Lagrangian framework, 
as the fit results can be affected in an unpredictable way, some action has to be taken. The simplest solution 
is certainly to discard the faulty data samples; however, as suggested by Bbl . a down-weighting of the outlier 
contributions to the minimized might also be considered. This could be a way to reconcile the preservation of 
the fit information quality with the use of all available samples. 
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4.1 The Basic x^/Least Square Method 

Usually, performing a fit - global or not - requires to minimize a funetior0 relating the 
differences between the measurements (m = {m,, i = 1 , • • • n}) and the corresponding model 
(theoretical) expections (M{a) = {Mi{a),i = 1, • • • n}) weighted by the error covariance 
matrix V provided together with the data spectrum. Leaving aside for now possible global (ad¬ 
ditive or mutiplicative) systematic uncertainties, the error matrix V provided by experimental 
groups gathers the statistical and systematic errors and, thus, is not necessarily diagonal. The 
vector a denoting the unknown internal model parameter list, minimizing : 

X^ = [m- M{a)fV-^ [m - M{a)] (2) 

with respect to a allows to derive its optimum value oq. When several independent data samples 
are to be treated simultaneously, the minimized is a sum of terms like Eq. Q, one for each 
data sample. 

As reminded in Pl5l . if the model M{a) is linear in the parameter^ and if the error co- 
variance matrix is correct, the estimated parameter vector do has unbiased components and 
this estimator do has the smallest variance. As illustration, in the case of a straight line fit 
(M = q + px), V. Blobel f[^ produced the residual plots for the model parameters using 
several kinds of error distributions for the generated data points (each with the same standard 
deviation) and showed that these plots are always gaussian distributions, as expected from the 
Central Limit Theorem. Of course, the probability distribution is flat only if the error distribu¬ 
tions are gaussian, i.e. if the effective function is actually a real x^- 

When analyzing (a collection of) actual spectra obtained by various groups, nothing better 
can be done and the derived fit solution faithfully reflects the whole data information on which 
it relies : It corresponds, at worst, to the least square solution and, at best, to the minimum x^ 
solution, depending on the functional nature of the true experimental error distributions. 

4.2 Iterative Treatment of Global Scale Uncertainties 

In the Subsection just above we have briefly summarized the traditional method which ap¬ 
plies when the handled spectra are not significantly affected by (correlated) global uncertain¬ 
ties. These can be of either kinds : additive (offset error) or multiplicative (scale/normalization 
error). As no offset error issue is reported for the spectra we analyze within BHLS [lMl 1^ . 
we skip this case and let the interested readers refer to suitable references W2\ l45l . In 
contrast, multiplicative (global scale) uncertainties are reported for most experimental spectra; 
when they are non-negligible compared with the other (more standard) kinds of errors, they 
should be specifically accounted for within the global fit procedure. This is of special concern 
for the important e+e“ —)■ vr+vr” data samples collected in scan mode [l52l l53l [54l . and even 
more for those collected using the Initial State Radiation (ISR) mode by KLOE [|^ ITTl |38]| . 
BaBar [|39ll40l or BESSIII [|4Tl : furthermore, the normalization uncertainties reported for each 
of the ISR data samples have all a peculiar structure which deserves each a specific treatment 
- this is the subject of the next Subsection. 

''Which is a true x^ if the errors are gaussian. 

'^Actually, fitting is generally performed in the neighborhood of some given solution; this makes the linearity 
condition less constraining in practice. 
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A constant global scale uncertainty, as those affeeting the data samples from CMD2, SND 
or BESSIII, ean be written (3 = 1 + where A is a random variable with range on ] — 1, +oo[. 
As E{\) = 0 and E{\^) = with a « 1, the gaussian approximation for A is safe [I43ll4^ . 
A data sample subjeet to sueh a global seale uneertainty provides an individual eontribution to 
an effeetive global whieh should a priori be written : 

\2 

= [m — M{d) — \AY'V~^[m — M(a) — XA] H —- (3) 

where m, M, V and a earry the same definitions as in Subseetion 14 .1 1 while A and a have just 
been defined. As for A, even if intuitively one may prefer A = m, the ehoiee A = M{a) has 
been shown to drop out any biasing issumi [I4^l45ll70ll . 

Assuming that the unknown seale faetor A is solely of experimental origin - and, then, 
independent of the model parameters a - the solution to dx^ /9A = 0 provides its most probable 
value Ao ll34ll . After substitution, Eq. ([3]) beeomes : 

X^ = [m — M{d)]^W~^[rri — M{d)] with W = V + a‘^AA^ (4) 

whieh exhibits a modified error eovarianee matrix W and only depends on the (physies) model 
parameters. More preeisely, the single reeolleetion of the seale uneertainty A is the oeeurenee 
of its varianee in the modified eovarianee matrix W. 

However, Eq. dH) elearly points toward a diffieulty if the model is not numerieally known 
beforehand as the modified eovarianee matrix beeomes a-dependent when setting the unbias¬ 
ing ehoiee A = M. In this ease, the parameter error eovarianee matrix provided by the x^ 
minimization might be uneasy to interpret. 

The way out is to define iterative proeedures; this is allusively stated in 021, but explieitly 
eonsidered in [[44l as solution to the so-ealled ’’Peelle’s Pertinent Puzzle’uj [|43l . provided a 
good starting approximate solution is known beforehand; however, defining sueh a tool might 
be a delieate task if the underlying model is non-linear, as quite usual in partiele physies. Sueh 
a proeedure has already been followed and sueeessfully worked out in WT\ in order to derive 
through a minimization proeedure the parton density funetions from several measured speetra. 
When dealing with samples of form faetor and/or eross seetion data, other appropriate iterative 
methods should be defined. 

The starting step of the iteration implies ehoosing some initial value for A, say A = Aq. 
Without further information, the best approximation one ean ehoose is obviously Aq = m, the 
experimental speetrum itself. Quite interestingly, this turns out to start iterating with A = 0 
(cr = 0 in Eq. dH)), i.e. (3 = 1, a unit seale faetor; this makes the eonnexion with the iterative 
method followed in [|47l . 

Then, the minimization of the x^ in Eq. dH) with A = Aq = m h performed using the 
MINUIT proeedure [17111 whieh yields the (step # 0) solutior0 Mq via the fitted parameter veetor 
value Oq. The next step (# 1) eonsists in minimizing Eq. dll) using A = Mq = M(do) whieh 

'^This does not mean that the choice A = m necessarily leads to a significantly biased solution as shown below. 

'"^Peelle’s reference is no longer of common access, but its main content - which closely resembles the 
D’Agostini issue raised in ll42l - is reproduced in li44l . 

'^The analysis method in ll34l[35l actually stops there; the present analysis aims at going beyond. 
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is easily implemented in the proeedure and, at eonvergenee, minuit provides the step # 1 
solution M(ai). This stepwise proeedure0. is followed until some eonvergenee eriterium is 
met. As in eaeh minimization proeedure the eovarianee matrix is eonstant, the interpretation of 
the parameter error eovarianee matrix is eanonieal. 

The eonvergenee speed of the iterative proeedure eannot be guessed ab initio but may be 
expeeted fast, referring to the fit of the parton density funetions where the eonvergenee is 
essentially reaehed at the first iteration [|47l . This is eonfirmed by the Monte Carlo studies 
reported in Appendix lAl 

Nevertheless, one may infer that the number of iteration steps is smaller for a starting guess 
for A elose to the aetual model than for an arbitrary ehoiee; elearly, as the ehoiee A = m (the 
experimental speetrum) should be the elosest to the aetual model, one may think that it should 
minimize the number of interations needed to reaeh eonvergenee. Additionally, this ehoiee 
does not imply any a priori assumption on the parameter veetor to be fitted. 

Among the data samples one deals within the BHLS based global fit method, most have 
been eolleeted in sean mode, essentially at Novosibirsk, and earry a eonstant seale uneertainty 
merging several effeets. This is espeeially the ease for the e+e“ —)■ tt+tt" data samples eol¬ 
leeted by the CMD2 [l52l|53l and SND [[5^ deteetors; this also eovers the ease of the BESSIII 
data sample iHTI . 

In order to simplify and unify the notations in the following diseussion, it is suitable to 
perform the ehange of random variable X = ap. Then, the statistieal properties for A propagate 
to E{p) = 0 and E{p?) = 1 and, defining in addition B = a A, Eq. (|3]) above beeomes : 


= [m — M{d) — pBYV ^[m — M{d) — pB] + 

The eondition dx^/dp = 0 provides the most probable value for p : 

B'^V-Am- Mid)] 

^ ~ B^V-^B + 1 

and, substituting this into Eq. ([5]), one gets : 

= [m-M{d)fW-^[m-M{d)] with W = V + BB^ 


(5) 

( 6 ) 

(7) 


Stated otherwise, from the point of view of the physies model, the minimization proeedure 
keeps traek of the seale dependenee by a modified eovarianee matrix whieh, in turn, influenees 
the fit. A faithful graphieal eomparison of data and model - like the usual fit residual plots - 
should take into aeeount the fitted seale, as illustrated in ll^ for instanee. 


4.3 Global Scale Uncertainties Effects in ISR Experiments 

With the advent of the faetory in Eraseati, of the J/ip faetory in Beijing and of the B fae- 
tories at SEAC and KEK, the possibility opened to get large data samples for the various e+e“ 
annihilation ehannels in the region of interest of the BHES model, namely, from the thresholds 
to the <p meson mass energy region (y^ < 1.05 GeV). The produetion meehanism involved is 

'®Each such step is defined as a full (MINUIT) minimization procedure where the covariance matrix is un¬ 
changed until convergence is reached. 
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the emission of a hard photon in the initial state iTT^ . the so-called Initial State Radiation (ISR) 
phenomenon. This ISR production mode has been used to collect high statistics data samples 
for the vr+vr” channel covering the low energies by the KLOE |l36l [33 EH, BaBar 

[[391I401 and BESSIII IHTII Collaborations. 

However, it is a common feature of the KEOE and BaBar (ISR) data samples to carry 
non-trivial error structures. Beside a non-diagonal statistical error covariance matrix (V), they 
exhibit a large number of (statistically independent) bin-to-bin correlated uncertainties, most 
of these being additionally s-dependent. As far as we know, this seems to be a premiere 
in particle physics and how this is dealt with inside minimization procedures deserves to be 
clarified and explicitely stated (see also ll35l l. 

Eet us consider a given experimental data sample E, a spectrum m function of s, for which 
the (given) statistical error covariance matrix is V ; the information provided for the bin-to-bin 
correlated uncertainties defines several independent scale uncertainties Aq (a = 1 , • • -nscaie) 
and should be understood as follows : each of the scale uncertainty Aq, is a random variable of 
zero mean and carrying a s-dependent standard deviation (Tq,(s) as tabulated by each experi¬ 
ment. It is clearer to make the change of (random) variables Aq, = (Ta(s)/iQ, (a = 1, • • • nscale) 
and assume that all the random variables fia fulfill E{jjia) = 0 and E^jjLajjLjs) = Sajs- 

Then, the other notations being identical to those previously defined, the iu Eq. ([5]) 
generalizes to : 

= [m- M{d) - - M(a) - (8) 


where implicit sum over repeated greek indices is understood. One has defined Ba{s) = 
crQ(s)y4(s), A being the s-dependent vector already defined. A is iteratively redefined as em¬ 
phasized in the previous Subsection. Using the minimum conditions dx^/d^a = 0 and the 
independence conditions of the various sources of scale uncertainty dfia/dfip = the most 
probable values for the /x„’s can be derived ll35ll . A recursion can be defined and allows to 
derive0 from Eq. dH : 


r X^ = [m — M{d)]^W ^[m —M(a)] 

(9) 

[ Wij = Vij + BiBj = Vij -f cr„(si)cr„(sj)] AiAj , (V[i, j]) 

in close correspondence with Eq. ©• 

A specific feature of Eq. (|9|) deserves to be noted. As each experimental group reports 
separately on each identified independent source of (scale) uncertainty, these should indeed be 
fitted separately as stated just above to go from Eq. dH to Eq. (jH). More precisely, for the 
experiment E, we are not using the quadratic sum (ct£;(s))^ = Z]a[c’‘a('5)]^ for its partial 
which would have given aE{si)aE{sj)AiAj inside the full error covariance matrix instead of 
what is shown in Eq. dH)- Stated otherwise, the various sources of normalization uncertainties 
are not summed in quadrature but really treated as statistically independent. 


4.4 Numerical Tests of the Global Fit Iterative Method 

As stated in the header of the present Section, if the physics correlations predicted by the 
Effective Eagrangian (here BHES) are fulfilled by the data, the estimate of the model parame- 

For clarity, defining Z = A or B, Zk denotes the quantity Z{sk) for short. 
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ters and the parameter error eovarianee matrix are legitimate tools to serve evaluating related 
physies quantities. 

As in the previous studies relying on the HLS model, at the early stages [|^ |33l or more 
reeently ||34l[35l|56l|571, the method is to minimize a global expression taking into aeeount 
the largest possible number of data samples and using appropriately all information provided 
by the experimentalists eoneeming all kinds of uneertainties whieh affeet their data samples. 
The aim of Subseetions 14.11 - 14.31 was to detail how the x^ pieee assoeiated with eaeh data 
sample should be eonstrueted, depending on its reported error strueture. 

In eontrast with previous referenees (ineluding ours), the fit proeedure will be adapted in 
the present study in order to examine and eure possible biases produeed by having stopped the 
fit proeedure at the A = m step instead of iterating further on as suggested in [l42l . explieitly 
proposed in [|44l and performed in ll47l . 

In order to eheek whether estimates based on global fit results ean be trusted as, for instanee, 
the muon HVP eentral value and its uneertainty derived from the fit information returned by 
MINUIT, some additional eheeks on the fitting method and its iterative aspeet deserve to be 
performed, at least to eontrol that, indeed : 

• The fit parameter residuals A* = af** — are unbiased gaussians, 

• The parameter pulls are eentered gaussians of unit standard deviations. 

One should also eheek that the fit probabilities distributions are uniformly distributed on [0,1] 
when the measurements are indeed true unbiased gaussian distributions. 

This eondition list ean be supplemented with some examination of the effeets due to non¬ 
linear dependenees upon the parameters to be fitted. 

However, eheeking this list of properties obviously implies that the true parameter values 
are known, that the measurements are indeed sampled on truely eentered gaussian distributions 
and that their errors are indeed the true standard deviations of the measured speetrum. Stated 
otherwise, this exereise goes beyond using aetual measured experimental data samples as, then, 
truth is unknown : The global fit method -as any other method - should be evaluated using 
data samples generated by Monte Carlo teehniques; in this ease, the true parameter values and 
their uneertainties are known at the sample generation level and ean reliably be eompared to 
the fit results. The detailed study is transferred to Appendix the most involved results are 
summarized right now : 

• The effeets of non-linear parameter dependenee within models used to fit data speetra 
(see Subseetion lA.2.11) are likely to be marginal for the kind of experimental distributions 
we are dealing with. This should be related with the loeal minimum finding strueture of 
the algorithms gathered within the MINUIT paekage. 

• When seale uneertainties dominate the sets of speetra globally submitted to fit, using0 
A = tue gives a solution whieh ean exhibit strong biases, but this solution is the start of 
an iterative proeedure whieh leads rapidly to the unbiased solution to the minimization 
problem. The biases oeeuring at start of the proeedure ean be very large, but they are 

ruE being the experimental spectrum in the expression for the (see Eqs. (I3]l or @). 


11 







observed to practically vanish already at the first iteration step (the previously called Mi 
solution). 

• When performing a global fit of some data samples dominated by global scale uncer¬ 
tainties together with others where the statistical errors (e.g. affecting randomly each 
bin) dominate, the iterative method obviously works as well as just stated. In this case, 
however, the presence of some samples free from scale errors exhibits an unexpected 
pattern : Even if the data samples free from scale uncertainties are affected by enlarged 
statistical errors, they strongly reduce the biases generated by the A = rriE choice. Stated 
otherwise, the effects of data samples where the normalization errors are dominated by 
the (random) statistical errors is to favor the smearing out of the biases in the parameter 
value estimations. 

The properties just listed concerning the unbiasing of the fit parameters extend to the esti¬ 
mates of physics quantitites derived from using the fit result information (parameter values and 
error covariance matrix). Additionally, as the parameter pulls are observed as centered gaus- 
sians of unit standard deviation, the calculated uncertainties relying on Monte Carlo sampling 
of the fit parameter distributions should also be reliable. This is of special relevance for the 
evaluation of the various contributions to the muon LO-HVP discussed in Section [3l 

The last item in the list just above has important consequences while working with real (and 
so, not really perfect) experimental data. However, even if the fraction of data samples free 
from - or marginally affected by - scale uncertainties may look large enough, it is nevertheless 
cautious to ascertain that the fit solution is indeed unbiased by performing one or two additional 
iterations. Indeed, the studies reported in Appendixl^tell that, anyway, the iterated fit solutions 
are always unbiased. 

Therefore, one may conclude from this Section and from the simulation studies reported in 
Appendix 1^ that global fit methods can indeed be trusted. The single proviso is that iterating 
the fit procedure as explained above is mandatory or, at least, cautious. 

The issue is now to examine how the results given in [[3^ and [|35l are modified when iter¬ 
ating beyond the approximation Ae = mg for all data samples significantly affected by scale 
uncertainties, constant (as, mostly, the spectra reported in [l52l l53l [54l ) or s-dependent (as all 
the ISR spectra reported in [[361 [371 EBl |39l). Observing the stabilizing effect of the data sam¬ 
ples dominated by statistical errors (like the and yr] final states) is also methodologically 
relevant. 

5 BHLS Global Fit Method : Present Status and Corrigen¬ 
dum 

As stated several times above, the Effective Eagrangian Model we use is the broken HES 
(BHES) model developped in Il34l . In this Reference, the BHES model is also applied to all 
data samples collected in scan mode, by the various Collaborations which have run on the 
successive Novosibirsk e+e“ colliders. These e+e“ annihilation samples cover the tt+tt", 7r°7, 
77 ,7r+7r“7r°, K~^K~, final states and have been discussed in detail in several previous 

studies [[3^ [331 [Ml; for the sake of conciseness, we will not repeat this exercise here. As 
the BHES model also covers the r decays from the early stages of its formulation [[27l . the 
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previous studies inelude the dipion spectra collected in the —)■ -K^Ti^Vr decay mode by 

ALEPH [[49117^ . Belle lISTl and CLEO Il50l . Also included within the BHES fit procedure are 
some light meson decay partial widths not connected with the annihilation channels already 
listed, like — ?• K*^ K^'y, rj' —>■ cey or 0 —)■ 

A second step has been to extend the study in [|34l to treat the high statistics ISR data 
samples for e+e“ —>■ tt+tt"; this has been the purpose of the study in [|35ll where the KEOE08 
[f3^ and KLOEIO ifJTIl data samples collected by the KLOE Collaboration and the data sample 
produced by BaBar ll39l have been examined. Since then, two new samples have been produced 
by the KEOE (KEOE12 BSl l and BESSIII iHTll Collaboration^ Except otherwise stated, all 
the fit results presented in this paper have been obtained using the Configuration B [[34ll {i.e. 
dropping out from the fit procedure the three pion data samples collected in the 0 mass region). 

The studies covered by [l34l [35l [561 ISTll rely on minimizing a global function summing 
up partial each associated with a given data sample. Eor each of the ~ 40 50 data 

samples, the partial was (canonically) constructed following the rules detailed in Section |4l 
However, as the fit was not iterated in the studies [|34l[35ll . it is worth checking to which extent 
the value of the muon HVP derived herefrom is changed by the iteration procedure. 

Eor the present study, a few coding bug fixes have been performed and a piece missing in 
the expression for the e+e" —)■ and —)■ yy cross sections has been included. So, 

when different, the results in the present paper supersede those in Il34l[35ll . 

As for the missing piece just mentioned : In the amplitudes y* —)■ yPo (Eq. (65) in Odll ) 
and the cross section formulae —)■ yPo (Eq. ( 68 ) in Odil ). the non-resonant piece should 

be modified as follows : 

(1 - Ci)Lp„ ^ (l - Lp„ . (10) 

This implies that the single process which depends separately on the EKTUY lITSl parameters 
C 3 and C 4 is the —)• 7 r+ 7 r“ 7 r° annihilation. In this case both C 3 + C 4 and^ — C 4 combinations 

come in, while all others quantities only involve the C 3 + C 4 combinatiorcj. We apologize for 
the inconvenience. 


6 BHLS Global Fit Method : Iterating with NSK Data Only 

In this Section, we report on global fits using the data reminded in the preceding Section 
and discussed in [f34l : as for the pion form factor data, we focus for the present exercise on 
using only the most recent scan data collected by CMD2 and SND [|74l[5^[53l|54l, excluding 
the older data samples from OEYA and CMD ITT^ . 

The CMD2 data samples are reported to carry constant bin-to-bin correlated uncertainties 
of 0.6% ( Il74l ). 0.8% ( [|5^ ) and 0.7% ( ll53l f. while SND reports a 1.3% constant scale uncer¬ 
tainty [|54l - except for their first 2 data points where it is 3.2%. Eor these data samples, the 

'^The KLOE12 and KLOE08 data samples are tightly correlated; actually, they mostly differ by their respec¬ 
tive normalization procedures. Comparing their respective behaviors within our global treatment is, therefore, 
interesting. 

^®The studies ll34l (351 have been performed hxing C3 = C4. The BHLS fit recovers a good ht quality by 
modifying the value for ci — C 2 as will be seen below. 
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XVN 

A = m 

Iteration Method 

A = M 

varying 


IMI 

A = Mo 

A = Ml 

Astart = Ml 

Astart ~ Mx 

Decays 

8.16/10 

8 .01/10 

8.03/10 

8 .01/10 

8 .02/10 

New Timelike tt+tt" 

121.54/127 

121.75/127 

121.75/127 

121.74/127 

121.75/127 

7r°7 

63.84/86 

63.98/86 

63.96/86 

63.98/86 

63.96/86 

VI 

120.87/182 

120.84/182 

120.84/182 

120.84/182 

120.83/182 

A 

+ 

1 

o 

101.82/99 

102.49/99 

102.43/99 

102.49/99 

102.43/99 

K+K- 

29.87/36 

29.77/36 

29.78/36 

29.78/36 

29.78/36 

K^lf 

119.21/119 

119.21/119 

119.18/119 

119.20/119 

119.19/119 

ALEPH 

19.67/37 

19.73/37 

19.71/37 

19.72/37 

19.70/37 

Belle 

28.24/19 

28.27/19 

28.29/19 

28.27/19 

28.29/19 

CLEG 

34.96/26 

34.82/29 

34.82/29 

34.84/29 

34.84/29 

X^/dof 

648.16/719 

648.85/719 

648.78/719 

648.85/719 

648.78/719 

Global Eit Probability 

97.2% 

97.1% 

97.1% 

97.1% 

97.1 % 


Table 1: Global fit results derived by using only the data from lf52l[53ll54ll for the e’''e —)■ 
TT+TT” annihilation. See the diseussion and eomments in SeetionO 
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partial are essentially given by expressions like Eq. dH). For the other data samples, we 
performed as in [l34l . 

The first data eolumn in Table [U displays the results of the fit performed by setting A = m 
in the assoeiated with eaeh experimental data speetrum generieally named m. The form 
faetor returned by this {A = m) global fit is named Mq and is used to perform the first iterated 
(^4 = Mq) global fit; the results of this fit are shown in the data eolumn #2; this iteration 
#1 global fit returns the solution named Mi. The iterated #2 fit is then performed by setting 
A = Ml in the expressions of the pion form faetor data samples, leading to another (M 2 ) 
solution; the fit results are displayed in the third data eolumn in Table [IJ 

One elearly observes a quite tiny ehange in the first iteration : 0.2 unit in the value of the 
7 r'''7r“ data samples; also the global ehanges by only 0.7 unit. When going from the first to 
the seeond iteration, the ehanges are almost invisible. This eorresponds for experimental data 
to the effeet reported in Subseetion lA.2.3l for our Monte Carlo data. As derived quantity, let us 
report on the leading order (LO) eontribution a^(7r7r) derived by integrating Eq. ([T]) between 
0.63 GeV/e and 0.958 GeV/e; using obvious notations, the previously reported fits yield : 

( A = m : a^(7r7r, [0.63,0.958]) = 358.95 ± 1.63 

I A = Mo : a^ivrvr, [0.63,0.958]) = 360.00 ± 1.78 (11) 

[ A = Ml : a^ivrTr, [0.63, 0.958]) = 359.99 ± 1.79 

in units of 10“^°. So, one observes a tiny effeet while iterating onee (0.3% for the eentral 
value) and no effeet when iterating twiee. In the present ease, where the former data from 
If75l have been dropped out from the fit, the ’’experimental” estimate is a^(7r7r, [0.63, 0.958]) = 
361.26 ± 2.66 (see Table 7 in [[34l f. 

Another way to aeeount for the seale uneertainty is to set A = M{a) (whieh depends on 
the parameters under fit) and perform the fit. A starting value for A must be ehosen (denoted 
Astart) but its value ehanges at eaeh step of the minimization proeedure. In this ease, the fit 
eonvergenee time is mueh larger than previously but the results are almost identieal to those 
already obtained by iterating. The last 2 eolumns in Table [T] display the fit results starting 
with Astart = Ml and also those starting from the fit solution derived herefrom (denoted MA- 
As for a^(7r7r, [0.63, 0.958]), the values derived in these last fits numerieally eoineide with the 
iterated eases displayed above. 

Therefore, one may indeed eonelude, as ean be inferred from the Monte Carlo studies 
reported in Appendix that the HVP value reaehed without iterating is very elose to the HVP 
derived from the onee iterated solution. One also observes, as expeeted, that iterating only 
onee already leads to the final result; indeed, from iteration #1 to iteration #2, the ehanges for 
a^(7r7r) are at the level of a few 10“^^. 

As for the fit quality refleeted by the x^ values at minimum and the eorresponding fit prob¬ 
abilities, the last line in Table [U indieates that, whatever is the way one treats the veetor A, 
they are all alike. This, onee more, eorresponds to expeetations, as ean be eheeked with the 
diseussion in Subseetion IA.2.31 and expeeially the properties of Figure [8l Nevertheless, it is 
useful to eheek that the twiee iterated solution does not modify the result derived from the onee 
iterated solution in a signifieant way. 
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7 BHLS Global Fit Method : Iterating Scan and ISR Data 

It remains to introduce the other tt+tt” data samples collected at colliders using the 
ISR mechanism. Reference lf35ll has already done this work with the data samples then avail¬ 
able using the method described in Subsection 14.31 without, however, iterating the procedure. 
The conclusion reached was that the KLOE08 [|^ and BaBar |[^ data samples have difficul¬ 
ties to accomodate - within the BHLS framework - the whole set of data samples covering the 
channels already reminded in Section [51 In contrast, the KLOEIO [[371 data sample was found 
to fit well the BHLS expectations. Complementing preliminary works |[56l[571 . we revisit here 
the issue with the two new data samples provided by KLOE (KLOE12) and BESSIII. 

7.1 The r+PDG Analysis 

In Ref. [[35l . it has been shown that the BHLS fitter can be run without explicitly using def¬ 
inite —)■ 7r’''7r” data samples besides the non vr+vr” channels. Indeed, on general grounds, 

one expects that some limited isospin breaking (IB) information specific of this annihilation 
channel can make the job together with the r dipion spectra. It has been shown that the par¬ 
tial widths r(ci;/0 —)■ tt+tt") and r(p° —)■ e+e“), together with the products {V = u, 4>) 
r(l/ 7r+7r“) X r(l/ e'^e~) represent an amount of information sufficient to reconstruct 
- within BHLS - the pion form factor in the channel. 



0.7 0.8 0,7 0.8 0.7 0,8 ^Q0y^ 


Ligure 1: The r-i-PDG prediction (red curve) of the pion form factor in annihilations in the 
p — oj interference region. The various superimposed data samples are not fitted; also displayed 
are the average distances of each of the —>■ tt+tt" data samples to the common r-i-PDG 

prediction. 

Before going on, it deserves noting that the decay information used to run the r-i-PDG 
method has been extracted from the Review of Particle Properties (RPP) [[bTl and that the above 
mentioned pieces of information are in no way influenced by the data collected by KLOE, 
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BaBar or BESSIII; actually, they are almost 100 % determined by the data samples from the 
CMD-2 and SND experiments. On the other hand, the r+PDG analysis is not influenced by 
the global scale issue whieh mostly motivates the present work. 

We have performed the r+PDG run using all annihilation data mentioned in the above 
Seetions (eonfiguration A [iMlO . The fit returns Xt/^t = 82.1/85 = 0.97. The best fit solution 
allows to reeonstruet the predicted invariant mass distribution of the pion form faetor in the 
TT+TT” annihilation; this predietion is expeeted valid over the whole BHLS range as 
shown by Figure 2 in Il35l . It is worth showing here the mass range from 0.70 to 0.85 GeV; 
Figure [U displays the r+PDG predietion on this range together with the available 7 r’'' 7 r“ data 
superimposed {and not fitted)', we have ealeulated the yfi distanee of eaeh sample over its full 
rangeEH. The average x^ per data point is indieated inside the eorresponding pannel. 

Figure [Uindieates that the average x^ distanees for the NSK (CMD-2 & SND), KFOFIO, 
KFOE12 and BESSIII samples are small enough to elaim a sueeess of the r+PDG method. 
One ean eonclude that they fulfill the eonsisteney issue diseussed in Seetion[^with the full set 
of data and ehannels covered by BHFS. One should note that the deseription of the BESSIII 
sample (whieh is not a fit) is as good as the fit published by the BESSIII Collaboration [|4T1l . 
For KEOE08 and BaBar, we reaeh the same eonelusion as in (SSI; nevertheless, one ean now 
eompare the behavior the twirls samples KLOE08 and KEOE12 : We have X^kloeos = 4.8 
while x'^KLOEu = 1-2 elearly refleeting a better understanding of the error eovarianee matrix, 
while the eentral values are almost unehanged, as elear from Figure [IJ 

Stated otherwise, the issue met with KEOE08 and BaBar is eonfirmed but the two new data 
samples published sinee If^ are both found is good eorrespondenee with expeetations. 

7.2 The Iterative Method : Global Fit Properties 

The issue is now to report on the behavior of the global fits performed using the iterated 
method when the tt+tt" ISR and sean data are eonsidered simultaneously; this eomplements 
the work already presented in Seetion when using the sean data only. Exeept otherwise 
stated, the r data samples are always ineluded into the fit proeedure. On the other hand, as the 
behavior of the global fit for data/ehannels other than tt+tt” does not differ sensitively from the 
information already displayed in Tabled! this will not be repeated. 

Table |2] displays our main results using the sean and ISR e+e“ —)■ tt+tt" annihilation data. 
They eorrespond to the iteration # 1 fit (denoted above A = Mq), however the previously ealled 
A = m or A = Mi solutions gives almost identieal fit quality result^ 

The first data line displays the global fit properties with the indieated e+e“ —)■ 7r’''7r“ data 
samples used eaeh in isolation within the global BHFS eontext, together with all other data 
samples eovering the rest of the eneompassed physies (see Seetion [5]). 

One observes that the average (partial) per data point is of the order 1 

or (mueh) better and the probability high when running with any of the KEOEIO, KEOE12, 

For BaBar, the computed referred to here is computed on its spectrum up to 1 GeV, but truncated from 
the drop-off region (0.76+ 0.80 GeV). 

^^They mostly differ by the normalization method used to reconsUuct the spectrum from the same collected 
data. 

As regard to the fit parameter values and uncertainties : The A = Mq and A = Mi solutions differ unsignifi- 
cantly; the A = m exhibits some small departure commented below. 
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Eit Configuration 

Iteration Method 


KEOE08 

KEOEIO 

KLOE12 

NSK 

BESSIII-III 

BaBar 

N + - 

(60) 

(75) 

(60) 

(127/[209]) 

(60) 

(270) 

Eits in Isolation 

1.64 

0.96 

1.02 

0.96[0.83] 

0.56 

1.25 

Global fit prob. 

59% 

97% 

97% 

97%[99%] 

99% 

40% 

Eit Combination 1 


1.02 

1.48 

1.18[0.96] 

0.56 

1.36(*) 

xl+^-lN^+^- & Gl. fit prob: 




1.21 & 22% 


Eit Combination 2 


1.00 

1.05 

1.11 [0.89] 

0.61 


xl+^-lN^+^- & Gl. fit prob: 



0.98 & 99% 



Eit Combination 3 


1.02 

1.05 

1.10[0.89] 



xl+^-/N^+^- & Gl. fit prob: 



1.06 & 97% 




Table 2: Global fit results as funetion of the —)• 7r’''7r“ data sample eontent. Each entry 

displays the value returned by the global fit. The data samples involved can be 

tracked from the column titles, the following line giving the corresponding data point numbers 
[iV^+^-] in the range up to 1 GeV. The value flagged by * has been obtained using a BaBar 
sample truncated from the energy region [0.76, 0.80] GeV (250 data points). 


NSB@ and BESSIII data samples; as in ll^ the picture is not as good for KLOE08 and BaBar. 

Performing a global BHES fit using the data samples from KEOEIO, KLOE12, BESSIII, 
NSK and BaBar (amputatecj^ from the energy region [0.76, 0.80] GeV) leads to results given 
at the entry lines flagged by ”Eit Combination 1”; as the correlations between the KEOE08 and 
KEOE12 data samples are strong and their content not explicitely statecQ it is more cautious 
to avoid dealing with the KEOE08 and KEOE12 samples simultaneously. Despite the removal 
of the drop-off region in the BaBar tt+tt” spectrum, the global fit quality looks poorer. 

The results obtained when using the KEOEIO, KEOE12, NSK samples within the fit pro¬ 
cedure are displayed at the Entry ”Eit Combination 2” when BESSIII data are also included 
and ”Eit Combination 3” when they are not; the data and fit corresponding to the ”Eit Combi¬ 
nation 2” are shown in Eigure [2l Both Eit Combination 2 and Eit Combination 3 are clearly 
satisfactory. 

here denotes the collection of data samples from CMD2 117411521 l53l . SND 15^ (127 data points in total) 
as well as the former (82 data points) samples collected by OLYA and CMD ifTSl . The numbers in Table|2lgiven 
within square brackets include the contributions from these former samples. 

^^We remind that this removal is motivated by a possible mismatch in the energy calibration in the — w 
interference region between BaBar and the other tt+tt” data samples submitted to the same global framework. In 
contrast, when running with the tt+tt” BaBar sample in isolation, its full spectrum is considered. 

^®Some work in this field seems ongoing ca. 
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Vs(GeV) 


Figure 2: The pion form faetor data and fit eorresponding to the iteration # 1 BHLS global fit. 
The —)■ 7r’''7r“ data samples are those shown in the entry ’’Fit Combination 2” in Table [2l 

The inset in the top panel magnifies the peak region. The dowmost panels magnify the 

behavior in both distribution wings. See Section|7^for further eomments. 


19 














Therefore, this proves that the sean data from CMD2 and SND are eonsistent with the 
KLOEIO, KLOE12 and BESSIII data samples and that all these are fully eonsistent with the 
other data speetra introdueed in the global fit proeedure as indieated by the global fit probabil¬ 
ity. One should also remark that the systematie uneertainties provided for KEOE12 lead to a 
satisfaetory global fit, in eontrast with KEOE08, as already noted in the previous Subseetion. 

Exeept otherwise stated, the fit parameter values presented from now on are derived using 
the e+e“ Tr~^Ti~ data samples eorresponding to the ”Eit Combination 2” (see Table |2l); the 
fit results are those derived after the first iteration and they do not differ signifieantly from the 
eorresponding results at iteration # 2. The fit quality for the non-7r’''7r” data samples are almost 
undistinguishable from the numbers already given in the seeond data eolumn from TablefH they 
are not repeated for the sake of brevity. 

7.3 The Iterative Method : Updating The Model Parameter Values 

Beside improving the fits by mean of the iterative method, the present work aeeounts for an 
error and a eouple of bugs affeeting our Moreover, the present work ineludes the new 

KEOE12 data sample within the fit proeedure; this is not harmless as KEOE12 eonstrains the 
fit eonditions more severely than the KEOEIO sample. Therefore, the present results update 
and supersede the eorresponding ones previously given in [l34ll^ . 

7.3.1 The HLS-FKTUY Parameters 

The non-anomalous HES Eagrangian (broken or not) ean be written : 


^HLS — A + O-HLS^V (12) 

The unbroken expression for Chls can be found in lfT3l and its broken expression (BHES) is 
given in lIMll . The eovariant derivative whieh allows to eonstruet both pieees of Chls intro- 
duees the fundamental parameter g, known as universal veetor eoupling. The eoeffieient uhls 
is a speeifie feature of the HES model, expeeted elose to 2 in standard VMD approaehes; how¬ 
ever, phenomenology rather favors uhls — 2.5 sinee the early applieations of the HES model 
to pion form faetor studies [177117^12^ . 

On the other hand, the anomalous (EKTUY) seetor [fTSll of the HES model [[T3l eonsists 
of 5 pieees (see also Appendix D in ll34ll l. eaeh weighted by a speeifie numerieal parameter 
not fixed by the theory. Using eommon notations IfTSl l34l and faetoring out, for eonvenienee, 
the weighting faetors, the EKTUY Eagrangian eolleeting all the anomalous eouplings ean be 

writteio : 

CpKTUY = C 2 ,CyvP + {Ci — C‘i)CAVP + {^~C/^)CAAP + {Ci — C 2 — C^)CvPPp-\-{Ci — C 2 -\-Ci)CAPPP 

(13) 

Actually, the erratum involved in Eq. Coll comes from having missed the contribution of the (c 4 — C 3 ) 
term displayed in Eq. CU which actually turned out to impose C4 = C3. As already stated, after correction, 
all the anomalous decay couplings and the amplitudes for e+e“ — (t'‘°/p)7 anihilations only depend on the 
combination (04 -f C3)/2 and the single place where the difference (04 — C3) occurs is the e+e" —>■ 
annihilation amplitude. In Il34l[35]| where (04 — C3) was absent, its physical effect was absorbed by (ci — C2) to 
recover good ht qualities; so (04 — C3) and (ci — C2) should carry an important correlation. 
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where P and V indieate the basic pseudoscalar and vector meson nonets and A the electromag¬ 
netic field. As Chls, ^fktuy depends on the universal vector coupling g. 

At iteration # 1, the global BHLS fit returns : 


(C4 + C3) 
“ 2 

= 0.956 + 0.004 

_ (C4 - C3) 
~ 2 

= -0.166 + 0.021 

Cl - C2 

= 0.915 + 0.052 

9 

= 5.507 + 0.001 

0-HLS 

= 2.479 + 0.001 


with correlation coefficients never larger than the percent level, except for < 5 g danis >= 
—0.30 and < 5 [ci — C2] (5[(c4 — c^)/ 2 ] >= +0.86. The sign of the {g, auLs) correlation term 
is easy to understand as the vector meson coupling to a pion or kaon pair rather depends on 
the product g' = auLsd- The large value of the ([ci — C2], [04 — C3]) correlation is also not 
surprising (see footnote [27]). The numerical values for g and anis are in the usual ball park 
and do not call for more comments than in lf34l[^ . 

Our value for c+ agrees with the estimates derived in [[T3ll from the 7r°77*) form factor 
(c+ = 1.06 ± 0.13) and from the u —)■ 7r°7 partial width (c+ = 0.99 ± 0.16) with a much 
smaller uncertainty due to the large amount of data influencing the (global) fit. After the bug 
fixing, c_ is found small but non-zero with a large significance and (ci — C2) becomes closer to 
1. Using the full 25 x 25 parameter error covariance matrix returned by the global fit, we have 
computed separately C4 and C3 by a Monte-Carlo sampling. This gives C3 = 1.124 ± 0.022 and 
C 4 = 0.789 + 0.021. 

Among the numbers displayed in Eq. (IT4l) . some are appealing : The nearness to 1 of 
the fitted ci — C 2 and c+ parameters, their customary guessed value ifTSll . should be noted and 
deserves confirmation with more precise data on the anomalous annihilations and light meson 
radiative decays than those presently available. 

7.3.2 The Iterative Method : Pseudoscalar Meson Mixing and Decay Parameters 

The BHLS symmetry breaking of the Lagrangian piece Ca leads to pseudoscalar physical 
fields constructed as linear combinations of their bare partners. The mechanism involved is the 
BKY mechanism extended so as to account for both Isospin and SU(3) symmetry breakings 
[|34l ; it can be complemented by the pseudoscalar nonet symmetry breaking scheme generated 
by the t’Hooft determinant terms [1791 . The main effect of these determinant terms is to provide 
the bare Lagrangian with a correction to the PS singlet kinetic energy term governed by a 
parameter A expected small (see Eq. (7) in [l34l f. 

The BHLS model connects to (Extended) ChPT [I25ll24l . especially its two angle 6*0 and 0% 
mixing scheme; in particular, it relates these angles to the singlet-octet mixing angle tradition¬ 
ally denoted 9p, together with the BKY breaking parameters 2 ;^, and to A [1341 . 

The upper part of Table [3] displays in its first data column our fit results in the general case. 
The fit value for 6^ is in good agreement with other expectations [l24l as well as that for 6q. The 
smallness of this has led us to impose 6^0 = 0 within fits which leads to the results shown in 
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General Fit 

Constrained Fit 

^0 

2.77° ± 0.41° 

0 

Os 

-25.95° ±0.35° 

-25.52° ±0.20° 

Op 

-15.29° ±0.32° 

-13.96° ±0.16° 

A 

(2.91 ±3.35) 10-2 

(1.86 ±1.17) 10-2 

^0 

(4.12 ±0.33) 10-2 

(4.00 ±0.33) 10-2 

e{rj) 

(5.85 ±0.48) 10-2 

(5.57 ±0.47) 10-2 

e'{7i') 

(1.46 ±0.13) 10-2 

(1.36 ±0.12) 10-2 

xVNdof 

Probability 

887.5/994 

99.3% 

892.5/995 

99.1 % 


Table 3: Some parameter values derived when leaving free Op and A (first data eolumn) or 
when relating them by imposing 6^0 = 0 to the fit (second data column). 


the second data column. The value for A undergoes a severe correction compared with [[34l[35]l 
and, presently, because of its large uncertainty, could be neglected without any real degradation 
in fit qualities. 

BHLS also allows for some additional contribution to the — rj — rj' mixing based on some 

possible aspects of Isospin breaking not already accounted for by the extended BKY scheme 

developped in [|34l . This turns out to redefine the physical (observable) fields (right-hand side) 

in terms of the (BHLS) renormalized (left-hand side) fields by [[^ : 

0 / / 

7 rR = vr - e 7]-e 7] 

< 7]\ = cos9p{7]-\-e tt'^) + sm.9p{7]'-\-e' tt'^) (15) 

7]p = — sin 9p{rj + e 7r°) + cos 9p{rj' + e' 7r°) 

Inspired by ifSOl . one can lessen the number of free parameters by stating : 


e = 


£ = 


CQ up- 


Vzcos 9p + sin 9p 

^ „ \/2 cos + sin 6*p 

-2eo sm 9p—^ - 

v2 cos 9p — sin 9p 


(16) 


and fit cq. Then, using the fit results (parameter central values and error covariance matrix), 
one can reconstruct the value for e and e'. The updated values are given in Table[3]still indicate 
a 7 r° — ?7 mixing much larger than the 7 r° — rj' mixing (a factor of 4). 
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Before closing this Subsection, we mention that the Monte Carlo sampling method allows 
to reconstruct the decay constant ratio fx/fir = 1-265 ± 0.009 which becomes fx/fir = 
1.295 dz 0.002 when constraining the fit with 9 q = 0. 


8 The Muon LO-HVP : Evaluations From Iterated Fits 

The main aim of the present study is to produce improved estimates of the muon LO-HVP 
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h-P-! 

T (A+B+C)+ PDG 




356.55±1.9g 

8 


: ■ : □ 

BaBar Amp. (0,28-1.00) 



i 

361.81 ±0.35(exp;365,2±2.7) 

7 

- 

i 

KLOE 08 



• ; * ; 

352,28±2.31(exp;356.70±3.13) 

6 

- 

i 

KLOE 10 



— 

353.22±2.35 (exp:353.30±3.26) 

5 

- 

hOji ! 

KLOE 12 



4-: i 

353.82±0.94 (exp:354.38±2.88) 

4 


• b—□—1 

CMD2+SND 



: — 

360.00±1.78 (exp:361.26±2.66) 

3 

- 

o : : 

KLOE(KL0E10 +KL0E12) 



TSTi : 

354.65±0.81 (exp:353.91±2.16) 

2 


o: ! 

NSK(all)+KLOE(10/12) 



: 

355.17±0.75 (exp;356.67±1.69) 

1 

- 

i : □ 

NSK+KLOE+BabBar Amp. 




358.37±0.39 (exp:359.07±1.43) 
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10'° a^(7T7T,Vs [ 0.63 -0.958 1) 

Figure 3: Values for [0.63, 0.958]) in units of 10“^° derived from global fits using the 

indicated —)■ tt+tt” data samples or combinations; the r dipion spectra are always used. 

The full green circles are the results obtained from the A = mfii (no iteration) and the black 
empty squares are the results obtained from the A = Mq fit (first iteration). The values de¬ 
rived by integrating the experimental spectra are indicated by red stars. See Subsection [8T| for 
comments. 


Il34l l35ll by means of the iterated global fit method expected to cancel out possible biasing 
effects which could affect the A = m {i.e. non-iterated) solution. The validity of the iterated 
method is supported by the Monte Carlo study outlined in Appendix which clearly indicates 
that the iterated method cancels out possible biases and returns, correctly estimated, the fit 
parameter uncertainties. Therefore, building on the conclusions collected in Subsection 14.41 
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one can produce bias free evaluations of the muon LO-HVP. The effects of iterating^ from 
Mq to Ml - the solution derived using A = Mq within the fit procedure - will be especially 
emphasized. To be complete, this update also takes into account the new KLOE12 [l38l and 
BESSIII W TT+TT data samples - which happen to be very constraining - and also corrects 
for some bugs. Therefore the present numerical results supersede the corresponding ones in 

mm. 

8.1 Various Evaluations Of a^(7r7r, [0.63,0.958] GeV) 

The point at top of Eigure|3]is the so-called r+PDG fl^ value for [0.63, 0.958] GeV) 

derived by switching off the contributions of the various —)■ tt+tt” data samples from the 

minimized replacing them by decay information extracted from the Review of Particle 
Properties (RPP) [|6T]| as emphasized in Subsection 17.11 

In order to get the other points displayed in Eigure [H one always uses all the channels 
covered by BHES, including the r spectra from AEEPH, CEEO and Belle. As for the —)■ 

TT+TT” data samples, one uses each of the BaBar, KEOE08, KEOEIO and KEOE12 samples in 
isolation as indicated within the Eigure (see also Table[2]and Subsection l7.2l) . The point flagged 
by CMD2+SND is obtained from a fit to the so-called [[34l new timelike data from CMD2 and 
SND dTll [5^ |53l |54l, leaving aside the older data from OEYA and CMD collected in ifTSl 
(see Tableland Sectionabove). As for the BaBar spectrum, for reasons already stated, 
the fit is performed on the spectrum amputated from the drop off region (i/s G [0.76, 0.80] 
GeV). Einally, as the published BESSIII spectrum ends up at 0.9 GeV, one cannot produce an 
experimental value on the interval [0.63, 0.958] GeV. 

As a general statement, Eigure |3] clearly illustrates that the iterated (Mi) and the non- 
iterated (Mq) solutions provide quite similar fit estimates for a^(7r7r, [0.63, 0.958] GeV). One 
should nevertheless remark that the agreement between both fit solutions and the numerical 
integral of the experimental data is less satisfactory for the data samples which exhibit poor 
fit qualities within the global framework (KEOE08 and BaBar) than for the others (KEOEIO, 
KEOE12, CMD2+SND) as can be inferred from the ’’fit in isolation” properties displayed in 
Tabled Einally, the weighted averages of the experimental results for KEOEIO and KEOE12 
alone or together with all NSK data (the so-called new timelike data and the former samples 
[iTSll l are always well reproduced by the global fit and are supported by quite good probabilities 
(see Table O. 

Using the NSK+KEOE(10/12) sample configuration, the iterated BHES global fit gives a 
slightly smaller central value (by ~ 1.5 10“^°) while the uncertainty is improved by a factor 
~ 2. It is also worth pointing out the role of the r spectra within the BHES global fit framework. 
The following numbers illustrate how the constraints involved by the r spectra allow BHES to 
yield a more precise fit estimate for a^(7r7r, [0.63, 0.958] GeV). Comparing the direct integration 
result to the values derived from fits, one indeed gets at iteration # 1 : 

r Direct Integration : a^(7r7r, [0.63, 0.958]) = 356.67 ± 1.69 

i A = Mo(fitexcl.r) : a/,(7r7r, [0.63, 0.958]) = 355.07 ± 0.96 (17) 

[ A = Mo (fit inch r) : a^(7r7r, [0.63, 0.958]) = 355.17 ± 0.75 

is the solution to the fit performed under the approximation already named in short A = m {i.e. each of 
the various tt+tt” experimental spectra is used for its individual contribution to the global x^)- 


24 







in units of 10“^°. 

Finally, the downmost point in Figure [3] displays the result derived using all data samples 
(except for KLOE08 as there is not enough published information to account for its strong 
correlation with KLOE12); this estimate for a^(7r7r, [0.63, 0.958]) which benefits from a very 
small uncertainty has, however, a poor fit probability as clear from Tabled 

8.2 Contributions To The Muon LO-HVP Up To 1.05 GeV 


Channel 

A = m 

A = Mo 

Exp. Value 

7r’''7r“ 

495.06 ± 1.43 

494.59 ± 0.89 

492.98 ±3.38 

7r°7 

4.53 ± 0.04 

4.54 ± 0.04 

3.67 ±0.11 

m 

0.64 ±0.01 

0.64 ±0.01 

0.56 ±0.02 

+ 

1 

o 

40.83 ±0.57 

40.84 ±0.57 

43.54 ± 1.29 

KlKs 

11.56 ±0.08 

11.53 ±0.08 

12.21 ±0.33 

K+K- 

16.79 ±0.20 

16.90 ±0.20 

17.72 ±0.52 

Total 

569.41 ± 1.55 

569.04 ± 1.08 

570.68 ±3.67 


Table 4: The contributions to the muon EO-HVP from the various channels covered by BHES 
from their respective thresholds to 1.05 GeV in units of 10“^° at start and after iteration. The 
last column displays the direct numerical integration of the various spectra used within BHES. 
The 7r+7r“ data samples considered are those flagged by ’’Combination 2” in Table[2l 

The EO-HVP’s integrated from their respective thresholds up to 1.05 GeV are displayed 
in Table S the central value for a^(7r7r) includes final state radiation (ESR) effects. The first 
data column shows the results from the fit solution Mq derived from fitting with A = m\ the 
second data column displays the results corresponding to the solution Mi derived by fitting with 
A = Mq. These two data columns report on the fits performed using all annhilation channels 
encompassed by BHES and the r dipion spectra. Einally, the rightmost data column provides 
the direct numerical integration of the experimental spectra - actually only those feeding the 
BHES fit procedure, including the KEOEIO, KEOE12 and BESSIII data samples besides the 
scan data. 

As for the 7r’''7r“ channel, both fits - which include the r spectra - provide central values in 
agreement with each other and with the direct estimate within the quoted erroiE^. If the A = m 
solution were (inherently) exhibiting a bias, comparing the first two numbers in the first line of 

As for the central value of the experimental estimate which is the present concern, one can legitimately expect 
that it should be affected by some bias (a priori, of unknown magnitude) of the same nature than the A = m result. 
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TablelHindicates that this does not exceed ~ 0.5 x 10“^° - e.g. half a standard deviation. There¬ 
fore, real experimental data samples confirm the gain provided by a global fit procedure when 
samples with normalization errors small compared to their statistical accuracies are included; 
exploring this effect is the purpose of Subsection lA.2.3l in Appendix lAl 

One should also remark that the unbiasing iterative procedure lessens significantly the un¬ 
certainty on a^(7r+7r“) compared with the A = m solution and, over the whole range of validity 
of BHLS (up to 1.05 GeV), one ends up with a factor of ~ 3 reduction of the uncertainty com¬ 
pared to the direct numerical integration. The same kind of effect is reported in Il47l concerning 
the spread of the parton density functionS 

Therefore, relying on the iterative procedure, one observes that the global fit does not pro¬ 
duce significant shifts of the central values of the HVP contributions which could be attributed 
to the normalization (scale) uncertainties strongly affecting some data samples. Relying on 
the Monte Carlo studies outlined in Appendix this can be attributed to the large number of 
data samples where the statistical uncertainties dominate over the normalization uncertainty. 
Moreover, the uncertainty on the part of the LO-HVP derived from the BHLS fit (more than 
80% of the total LO-HVP) is very small and even marginal. 

8.3 The Muon g — 2 From BHLS Global Fit Procedure 

In order to evaluate the muon LO-HVP from the fit results derived by means of the BHLS 
global fit procedure, the numbers given in Table |4] should be supplied with several additional 
contributions which cannot be derived from within the BHLS framework but should be esti¬ 
mated by other means. This covers the channels opened below 1.05 GeV but remaining out¬ 
side the present BHLS scopeEi] and, more importantly, all hadronic contributions covering the 
non-perturbative QCD region above 1.05 GeV should be estimated via the direct integration 
method. 

Table [5] summarizes these additional contributions to be combined with the BHLS results 
to derive the muon LO-HVP; in this Table, one reminds the information available by end 
of 2011 and used in our previous [l34l [35ll . The data column flagged by ’’LO-HVP (2014)” 
is the update derived by taking into account the data samples more recently collected (and 
published up to the end of 2014); these are the e+e“ —)• 3(7r+7r“) data from CMD-3 (Si], 
the —)■ cutt® —)■ 7r°7r°7 from SND [[82l and several data samples collected by BaBar in 

the ISR models [[83] l84l [85] |86]|. These data samples highly increase the available statistics 
for the annihilation channels opened above 1.05 GeV and lead to significant improvements. 
One thus should note the important improvement these provide for the LO-HVP contribution 
from the [1.05, 2.0] GeV region : its uncertainty is reduced by 25 %, while its central value is 
almost unchanged. Despite this improvement, the energy region [1.05, 2.0] GeV still remains 

Indeed, roughly speaking, the experimental cross section <7exp{s) is related with the underlying theoretical cross 
section ath{s) by a relation of the form aexp{s) = (Jth{s) -f 5a{s) and the 5a{s) correction depends on the 
normalization uncertainties which just motivate the iterative method! Actually, this Sa(s) is exactly the scale 
dependent term in Eqs. (|5]l and ® . Obviously it cannot be estimated without some htting procedure. 

^°In particular, Figure 5 in this Reference, is quite informative about the variety of correction kinds revealed by 
unbiasing procedures. 

^*For instance the 4, 5 of 6 pion annihilation channels, or the cutt^ final state. 

^^These cover the pp, K'^K~, K^Ks, KlK' n'~, KsKsTT~^Tr~, KsKsK'^ annihilation final states. 
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Contribution from 

Energy Range 

LO-HVP (2014) 

LO-HVP (2011) 

missing channels 

threshold —> 1.05 

1.34(0.03)(0.11)[0.11] 

1.44(0.40)(0.40)[0.57] 

J/ij 


8.94(0.42)(0.41)[0.59] 

8.51(0.40)(0.38)[0.55] 

T 


0.11(0.00)(0.01)[0.01] 

0.10(0.00)(0.01)[0.01] 

hadronic 

(1.05, 2.00) 

60.45(0.21)(2.80)[2.80] 

60.76(0.22)(3.93)[3.94] 

hadronic 

(2.00,3.10) 

21.63(0.12)(0.92)[0.93] 

21.63(0.12)(0.92)[0.93] 

hadronic 

(3.10,3.60) 

3.77(0.03)(0.10)[0.10] 

3.77(0.03)(0.10)[0.10] 

hadronic 

(3.60, 5.20) 

7.50(0.04)(0.05)[0.06] 

7.64(0.04)(0.05)[0.06] 

pQCD 

(5.20, 9.46) 

6.27(0.00)(0.01)[0.01] 

6.19(0.00)(0.00)[0.00] 

hadronic 

(9.46, 13.00) 

1.28(0.01)(0.07)[0.07] 

1.28(0.01)(0.07)[0.07] 

pQCD 

(13.00,oo) 

1.53(0.00)(0.00)[0.00] 

1.53(0.00)(0.00)[0.00] 

Total 

1.05 —)■ cxD 

+ missing channels 

112.82 ± 3.01toi 

112.96 ±4.13toi 


Table 5: LO-HVP contributions to with FSR corrections included. The statistical and 

systematic errors are given within brackets; the total uncertainty is given within square brackets. 
Column ’’LO-HVP (2011)” displays the contributions estimated using only the data samples 
available in 2011; Column ’’LO-HVP (2014)” displays the corresponding values updated with 
the data samples published up to the end of 2014. 
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the dominant uncertainty on the muon LO-HVP and this strongly limits the effect of gaining 
further in precision on the part of the LO-HVP covered by BHLS. 

Deriving the full HVP value also requires to account for the higher order effets. This 
includes the next-to-leading order contribution (NLO) taken from If26l ([—9.97±0.09] x 10“^°) 
and the recently estimated next-to-next-to-leading order (NNLO) effects which happen to be 
non-negligible ([1.24 ± 0.01] x 10"^°) If87l . 

To compute the muon g — 2, one should also include the light-by-light (LBL) contribution 
(here taken from [[88]|), the QED contribution If89ll90ll and the electroweak contribution (EW) 
ifSTI . The next-to-leading order contribution to the EBE amplitude (NEO-EBE) has also been 
computed recently If^ but is clearly negligible ([0.3 ±0.2] x 10“^°). Altogether, the numerical 
values we use (see Table [b]) are rather consensual [l92l . 


10^° X 

Value 

scan only 

s (incl. r) 

scan © KEOE © BESSIII 

Direct Integration 

scan © KEOE © BESSIII 

EO-HVP 

683.26 ± 3.78 

681.86 ±3.20 

683.50 ±4.75 

HO (NEO) HVP 

-9.97 ±0.09 ^ 

NNEO HVP 

1.24 ±0.01 ffsa 

EBE 

10.5 ±2.6 UMl 

NEO-EBE 

0.3 ±0.2 Ml 

QED 

11 658 471.8851 ±0.0036 [[89ll90ll 

EW 

15.40 ± O.lOhad ± 0.03Higgs,top,3-loop ttSIl 

Total Theor. 

11 659 172.62 ±4.60 

11 659 171.22 ±4.13 

11 659 172.86 ± 5.42 

Exper. Aver. 

11 659 208.9 ±6.3 

Attf, 

36.28 ±7.80 

37.68 ± 7.53 

36.04 ±8.31 

Significance (na) 

4.65ct 

5.00a 

4.38a 


Table 6: The various contributions to Aa^ = is given in units of 10“^°. For the measured 

value we have adopted the value reported in the RPP which uses the updated value for A = g-fi/gp recom¬ 
mended by the CODATA group ll93l . By KLOE, one means that the KLOEIO and KLOE12 tt+tt” data samples 
are introduced in the BHLS fit procedure and in the directly integrated spectra. 

The first data column in Table reproduces (after our methodological update) the muon 
anomalous moment estimate coming from the corresponding BHES global fit where only the 
scan data for the 7r’''7r“ channel are considered while all ISR data are excluded. This supersedes 
the corresponding information in |[34ll . The sample combination preferred by the BHES global 
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fit gives the results displayed in the second data column; it exhibits a 4.9a significance for a 
non-zero Aa^ = The evaluation derived by direct integration of the spectra used 

within the global fits are given in the third data column. The new data, as a whole, increase the 
discrepancy for Aa^ which is always found above the 4 ct level; effects of additional and not 
still accounted for systematics will be examined in the next Subsection. 

Figure |4] displays the results for Aa^ derived using or not the r data and various combi¬ 
nations of the available tt+tt” data samples introduced within the BHLS global fit procedure 
at first iteration. For comparison, one also displays in this Figure the evaluations produced by 
other authors and flagged by Dhea09 (H, DHMZIO [El, JSll (El and HLMNTll [j^ - 
corrected however for the recently calculated NNLO-HVP and NLO-LBL - contributions as 
included in Table [6l A priori, the Dhea09 estimate compares exactly to our evaluations us¬ 
ing scan data only; the other results are derived using, beside the NSK samples, the BaBar, 
KLOE08 and KLOEIO samples. These may compare to the last couple of lines in Figure 
m where the scan data are supplemented with the BaBar (not truncated), KLOE (10/12) and 
BESSIII samples. 

The following comments are in order here : 

• 1 / The difference between our estimates and those of other authors mainly concerns the 
estimated central value for Aa^. Also, our uncertainties are now reduced because of the 
global fit method, but also because of using much more data samples than other authors; 
this is clear by comparing the errors shown in EigurejUwith those given in [[35l . 

When using only the scan data, the shift one observes should reflect the biasing effect 
certainly present in the experimental data (see footnote # [29l) and corrected in our ap¬ 
proach by the iterated fit method. When the ISR 7r’''7r“ samples are also involved, the 
issue just reminded is amplified because the weight of samples with large overall scale 
uncertainties is much increasecj^. The effect of the BaBar data sample is no longer 
enough to balance the effect of the new data samples as clear by comparing the lines for 
”NSK-i-KEOE-i-BESSin” with the lines for ’’Global (ISR-i-scan)” which also include the 
(full) BaBar sample. Nevertheless, one should note the large difference of the correpond- 
ing probabilities. 

• 2 / When a comparison between a Aa^ estimate derived using the r data and the corre¬ 
sponding one excluding these is possible, ours exhibits the smallest difference (1.12 x 
10“^° for NSK-i-KEOE-i-BESSIII, —0.7 x 10“^° for the Global fit including all the 7r’''7r“ 
data samples). This is certainly due to the vector meson mixing which defines the BHES 
model. It is interesting to note that the IS 11 [l26l value, which is based on the 7 — 
mixing by loop transition^ is the closest to ours. 

• 3/ Relying on the global fit properties, the BHES model favors the ”NSK -i- KEOEIO 
-I- KEOE12 -i-BESSIII -i- r” as the largest consistent set of data sarnies. This leads to 
Aa^ = (37.55 ± 4.12) x 10“^° which exhibits a 5.a significancq^. Our estimate is 

^^All ISR data samples are strongly dominated by overall scale uncertainties, additionally s-dependent. 

^'^Within the BHLS model too, the 7 — mixing is mediated by loop effects. 

^^If using the data from 2011 in Table jSj as in our previous studies, this significance is ’’only” 4 . 8 ( 7 . This 
compares more directly to the results from other authors displayed in Figure (|4ji. The increased significance is a 
pure consequence of the recent improvements of the hadronic contribution from the [1.05, 2.0] GeV region. 
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expected to be free of biases generated by the overall scale uncertainties which dominate 
the ISR vr+TT data samples. 


r(A+B+C)+PDG 


NSK (CMD2+SND) 

KLOE 08 

KLOE 10 

KLOE 12 

BESS III 

BaBar (Trunc.) 

BaBar (Full) 


NSK (CMD2+SND)+r 
NSK 

DHea09 (e+e ) 


NSK+KLOE+BESS&r 
NSK+KLOE+BESS 
DHMZIO (e+e +t) 
DHMZIO (e+e ) 
HLMNTll(e+e-) 
JSll(e''‘e“ + r) 

Global (ISR scan&r) 
Global (ISR scan) 


BNL-E821(avrg) 


tti: : 


tea 

tea 




[35.30 

±4.58] 

[4.5 rr] 




Individual ttti 

' Data 

Sets + 

T - 


[35.97 

±4.63] 

[4.6 cr] 

[3/iVx. 

0.96] 

[99.5%] 

[38.78 

±5.16] 

[4.8 (t] 

[3/iVx. 

1.64] 

[58.9%] 

[39.21 

±5.15] 

[4.8 (t] 

[3/iVx. 

0.96] 

[96.6%] 

[38.33 

±4.33] 

[5.0 (t] 

[xVN,. 

1.02] 

[96.9%] 

[33.02 

±4.69] 

[4.2 7 


0.58] 

[99.9%] 

[29.15 

±4.07] 

[3.9 cr ] 


1.15] 

[73.8%] 

[27.40 

±4.03] 

[3.7 (t] 


1.25] 

[40.1%] 


0.96] 

0.97] 


— scan TTTT Data - 

I 7 V, [35.97 ±4.63] [4.6 ff] [xVA',. 

[37.94 ±4.95] [4.7 <t] [xVA^., 

[28.56 ± 5.8] [3.4 a] 

— scan +ISR tttt Data- 

[37.68 ±4.12] [5.0 cr] [xVlVx^ 0.90] 
[38.67 ±4.17] [5.1 cr] [x7W^„ 0.88] 
[17.96 ±5.4] [2.2 ff] 

[27.16 ±4.9] [3.3 ff] 

[24.56 ±4.9] [3.1 ff] 

[27.66 ± 6.0] [3.2 a] 

[37.02 ±4.03] [5.0 cr] 
tea [36.33 ± 4.03] [4.9 cr] 


[99.5%] 

[99.8%] 


[XVN.: 

W/N.: 


1.15] 

1.14] 


[99.1%] 

[99.7%] 


[18.5%] 

[24.3%] 


[0 ± 6.3] 


experiment 


-10 


40 


90 


140 


- “?) X10“ 

Figure 4: The deviation Aa^ = in units of 10“^°. The various have been derived 

from the global fit using the indicated e'''e“ —)■ 7r’''7r“ data samples and including/excluding 
the T dipion spectra as indicated. In red we display Aa^ corresponding to the iterated solu¬ 
tion and in green those corresponding to the A = m (non-iterated) solution. In blue results 
from other studies are given corrected by the recently evaluated next-to-next-to-leading order 
contribution IfSTl . See Section [83] for comments. 


8.4 Additional Systematics On The BHLS Estimate For The Muon g — 2 

A detailed study of additional systematics possibly affecting the BHLS evaluation of Aa^ 
has been already performed in |[35l . It concluded to an uncertainty of the LO-HVP central 
value for Aa^ = in the range [—1.3 0.60] x 10“^° coming from tt+tt” contribution 

in the 0 mass region, where BHLS is weakly constrained. An uncertainty coming from using 
the r spectra has also been considered; it was argued that the best motivated evaluation of this 
is the difference between fitting with the r spectra and without them in the most constrained 
configuration. Presently, this means that the BHLS preferred value (Aa^ = (38.58 ± 5.04) x 
10“^°) could be underestimated by ~ 0.9 x 10“^°. 

Another mean to detect systematics is to compare with the accurate ChPT predictions on 
the P-wave tt+tt” phase-shift [[94l and also with the available experimental data from the 
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Cem-Munich [[95ll and Fermilab ||9^ groups. These are shown in Figure [51 Included also are 
the predictions derived from the Roy Equations [|97l and from the phase of the pion form factor 
fit performed in ll26l (JSl 1). 

As for the BHLS predictions corresponding to using NSK+KLOE(10/12), we display in 
this Eigure the phase of the full amplitude and those corresponding to dropping out the isospin 
breaking (IB) effects due to the vector meson mixing^. The r spectra are included within the 
fit procedure. 




Vs (GeV) Vs (GeV) 


Eigure 5: P-wave tt+tt" phase-shift data and predictions from [l94ll (CGE) and ll26l (JSl 1) to¬ 
gether with the BHES phase-shift. The insets magnify the various behaviors close to threshold. 
See Subsection l8.4l for further explanations. 

The standard BHES phase shift predictions are displayed in the left-hand side panel of 
Eigure [51 One clearly observes a very good prediction of the phase-shift up to about 1.2 GeV, 
i.e. much beyond our fitting range (from threshold to 1.0 GeV for the tttt data). Indeed the 
Cern-Munich data are very well accounted for and the BHES predictions are in accord with 
the other predictions. The inset, however, exhibits a (minor) issue for the full amplitude phase, 
a small bump of about 1° close to threshold, absent from the IB amputated amplitude. This 
can be tracked back to a peculiarity of the broken HES model which does not split up the 
HK (Eagrangian) masses for the u and mesons and, consequently, the mixing angle a(s) 
does not exactly vanish at s = 0 (see Eigure 6 in [[321 ): in contrast the other angles fulfill 
/5(0) = 7 ( 0 ) = 0. Indeed, one has : 

fi(-s) 

where |[34l ei(s) is the difference of the charged and neutral kaon loops and n.,r 7 r(s) is the pion 
loop which both vanish at s = 0. This assumption has been checked with fits by imposing 

^^This is obtained by cancelling out the ’’angles” a(s), /3(s) and 7 ( 3 ) from the full amplitude expression. 
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= (1 + and choosing various fixed values for rj; the right-hand side panel 

in Figure [5] displays the phase shift for rj = 5% and, quite satisfactorily, its inset does not 
reveal a bump any longer. A non-zero (HK) mass difference rj cannot be generated 

by the breaking mechanisms already implemented within BHLS. However, a breaking of the 
nonet symmetry in the vector meson sector (VNSB) enables such an effect; this turns out to 
modify the customary vector field matrix - actually U(3) symmetric - within the covariant 
derivatives of the HLS model [IT3ll by a perturbation term proportional to the singlet vector 
field combination. The effect of VNSB has been derived from specific fit studies and indicates 
that Aa^ might have to be lessened by about 1.4 x 10“^°. 

Therefore, in total, the BHLS favored result can be expressed, in units of 10“^° as : 

Aa^ = 37.68 + ± 4.12*;, ± 6.3^.*, (19) 

where the three additional contributions play as shifts on the central value. Adding them up 
linearly, the maximum shift (—2.7 x 10“^°) may reduce the central value to 34.85 x 10“^° 
which has still a 4.6cr significance. The effect of these additional systematics is to reduce 
potentially by ~ 0.3cr all the significances displayed in Figure HI These are not due to overall 
scale uncertainties already accounted for by the iterative method; they might be reduced by 
new annihilation data samples covering the region up to 1.05 GeV in all the physics channels 
in the realm of BHLS. 


8.5 The HVP Slope At Origin In BHLS Fits 

In the lattice QCD approach of calculating extrapolation methods have been developed 
(see e.g. contributions to ll9^ ) to overcome difficulties to reach the physical point in the space 
of extrapolations. The low behavior of the euclidian electromagnetic current correlators 
on a lattice, which exhibits a discrete momentum spectrum, poses a particular challenge (see 
e.g. [|99l II0011 and References just below). The analysis of moments of the subtracted photon 
vacuum polarization function n(Q^) was particularly advocated in variants in Refs. UlOlll and 
[Il02l . Recent lattice calculations 01031110411105111061 have been utilizing moment analysis 
techniques for a more precise evaluation of The leading moment is given by the slope of 
the Adler function II1071I108I . the latter being given by : 


D{Q^) = 


R{s) 


-ds 


= —Q 

a dQ^ 


had I 




( 20 ) 


'Smin (S + 

where R{s) is the hadronic spectral functioi0 and Smin the smallest threshold energy squared 
(Smin = 'm'io within BHLS). Then, defining : 

Ris) 


Pi = 


■ ds 


( 21 ) 


the HVP slope at the origin is given by : 


d 


— Aahad(s) 
ds 


-s ^>+0 


37r Jsrr,.ir,. STT ^ 


( 22 ) 


^’i?(s) = cr(e+e — >■ hadr.)/a{e'^e — ^ yd'^ ) with cr(e+e —^ yd'y ) = 47ra^/3s by neglecting the 

electron mass. 
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The constant Pi can be directly estimated from data and partly from the BHLS fits. There- 
fore^ne can proceed as done above with our evaluations of and derive the results gath- 
erecoin Table |71 Here, one observes that the difference between the experimental and the HLS 
values for the HVP slope are at the percent level (a 2 aexp effect) and the uncertainty is scaled 
down by a factor of 10. However, to really feel the HLS improvement on the slope, one needs 
once more an improved hadronic spectral function at high energies. 


moment 

data direct 

HLS channels data 

HLS model 

HLS + non HLS 

Pi (GeV-2) 

11.83±0.08 

10.07± 0.05 

9.970± 0.016 

11.73 ±0.06 

as 

-0.92 ±0.01 

-0.78 ±0.01 

-0.772 ±0.001 

-0.907± 0.01 


Table 7 : The slope of the photon HVP at s = 0. 


A lattice estimate of the Adler function slope D'{0) has been presented in [I109II . The 
result is Pi = 5.8(5) GeV“^, and has been compared with Pi = 9.81(30) GeV“^, a result 
estimated using a phenomenological toy-model representation UllOII of the isovector spectral 
function. The lattice results too include the isovector part only and is missing higher energy 
contributions above 1 GeV. 

In the study I1102II . the authors provide numerical values from fits to Lattice data based on 
Fade approximants (PA). For this purpose, they parametrize the HVP as : 


n(Q2) = n(o) - 


N 

ao + Yl 

n=l 


Qjfi 

bn + 


(23) 


which thus leads to : 




(0) = 47ra 


dU 


( 0 ) 


—dvra 


N 

+ XI 

n=l 



(24) 


The parameters corresponding to the results they consider as optimal are given in their Table 3. 
Using their notations, their fitted parameter values leac j^ for instance, to (0.71 ± 0.15) x 10 ^ 
(PA solution [0,1]) or (0.75 ±0.30) x 10“^ (PA solution [1,1]). These compare reasonably well 
to the slope results reported in Table |7] just above, taking into account the proviso expressed 
above about lattice data. 


9 Concluding Remarks 

The present study was motivated by the question which gives its title to this paper. More 
precisely, the issue is whether the D’Agostini bias ||42]|46]| prevents to derive unbiased physical 
results from global fits to experimental spectra affected by dominant overall scale uncertain¬ 
ties 

^^The non-HLS part of Pi amounts to 1.76 ± 0.06 GeV“^. 

Assuming also the errors on the a’s and b’s parameters are not correlated. 

'*®We gratefully acknowledge G. Colangelo to have pointed out the issue for estimating the muon HVP using 
global ht methods. However, the bias issue is more general as will be argued shortly. 
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Actually, several issues are merged together. First, the effeetive global funetions to 
be used in the minimization proeedure should be appropriately defined. For the data samples 
where the statistieal errors dominate the overall seale uneertainties, the eonstruetion of the asso- 
eiated partial x^’s is quite standard. The real issue starts when the data samples are dominated 
by overall seale uneertainties. For eaeh of them, substantially, the eanonieal partial x^ has been 
reminded in Seetion[^and writes [I4^|45ll4^ : 

= [m — M{a) — \AY'V~^[m — M{a) — XA], 

leaving aside the so-ealled ’’penalty term” [|46ll proportional to A^. The (partial) x^ being 
appropriately defined, another issue is the ehoiee of the veetor A. 

In our former studies Il34ll35l . beside the ~ 40 data samples dominated by statistieal errors 
whieh follow the traditional treatment, the data samples eovering the —)■ 7r’''7r“ annihi¬ 

lation ehannel are all, sometime very strongly, dominated by overall seale uneertainties; this 
espeeially refers to the samples eolleeted by the KLOE and BaBar Collaborations using the 
ISR produetion mode. Here, for eaeh sample, we ehose for A the experimental speetrum itself; 
this ehoiee has been referred to as A = m all along the paper. The guess behind was that 
all seale uneertainties affeeting the different experimental speetra independently of eaeh other 
should smear out possible biases in the eentral values of the (eommon) theoretieal form faetor 
funetion parameters (SSI. 

It happens that the results one ean derive in this way from the BHLS global fit undergo 
very small biases (eompared to the errors derived from the fit proeedure); this is shown in the 
present stud)0. However, the guess just reminded was ineorreet and the aetual reason whieh 
explains the almost bias free results is following : As shown in the Monte Carlo study presented 
in the Appendix, there is no smearing out of biases if all the speetra submitted to fit undergo 
eomparable strong seale uneertainties; however, this study also shows that, if some of the fitted 
speetra are dominated by (random) statistieal errors rather than global seale uneertainties, the 
fit results ean be strongly unbiased. 

Nevertheless, a high level of unbiasing eannot be taken as granted as the real weight of the 
samples dominated by statistieal errors within the full global fit proeedure eannot be aseertained 
beforehand. Basieally, the ehoiee A = m potentially leads to biases of unknown magnitude; 
this has been shown by G. D’Agostini ll42ll with a simple example and more generally argued 
by V. Blobel Il4^ . These authors also showed that all biases vanish if, instead of A = m, one 
makes the ehoiee A = M, the ’’true” speetrum. But this is just not possible within eontexts like 
ours, where fits are performed just in order to derive the ’’true” speetrum from data. Fortunately, 
iterative methods allow to eireumvent this diffieulty by taking the path opened in WT\ in order 
to derive the parton density funetion from data and eorreet for biases. The iterative method 
we propose has been tested with the Monte Carlo study reported in the Appendix and shown 
to produee unbiased results with a quite fast eonvergenee speed; indeed, only one iteration is 
suffieient. 

So, our main eonelusion is indeed that global fit methods ineluding a fast iterative proeedure 
are expeeted to produee reliable pieees of information as, methodologieally, the eentral values 
are unbiased and the estimate for the uneertainties reliable; this espeeially applies to the part of 
the muon leading order HVP derived from annihilation eross seetions. 

^'which also corrects for some coding bugs affecting our previous studies. 
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Having shown that appropriate global fit methods should lead to results whieh ean be 
trusted, a related remark is worth being expressed. Iterative global fits allow to supply the 
BHLS Effeetive Lagrangian eross seetions with reliable and unbiased numerieal eentral values 
for the fit parameters and a good estimate of their error eovarianee matrix. Then, using these 
eross seetions and the fit information, Eq. ([T]) is expeeted to provide an unbiased estimate for 
a^( 7 r 7 r) as the ingredients are unbiased. 

On the other hand, when eomputing a^{ 7111 ) by direetly integrating a dipion speetrum in 
order to derive its so-ealled experimental value, one has to plug into Eq. ([T]) the experimentally 
measured eross seetion aexp.{s). However, as already noted in footnote #[291 or as ean be 
inferred from the eanonieal expression reminded just above, the experimental and model 
eross seetions are related by : 

^exp.(^) ^theorXs') T (5ct(s) 

where the best estimate of the seeond term write^ 5a{s) = XatheorXs). As obvious from Eq. 
db]), the best estimate of the seale faetor A equally depends on the measured speetrum and on 
the ’’true” speetrum, whieh ean be identified with its (iterated) fit solution. So, using again 
self-explanatory notations, Eq. ([B leads to : 

afj,{7m,exp.) = a^j_{Tni,theor.) + Sa^^rni) 

and thus a^( 7 r 7 r, exp.) looks intrinsieally biased for any sample subjeet to strong enough overall 
seale uneertainties. This issue is also refieeted by the residual plots whieh are improved when 
plotting the eorreeted residuals [m — (1 + A)M(a)] instead of the raw ones [m — M(a)], as ean 
be seen in Eigure 13 of ll35l : this allows to infer that is small but non-zero. It amounts 

to 5 a^( 7 r 7 r) ~ 2 x 10“^° in the ease ”NSK+KEOE10+KEOE12+BESSIII+r” favored by the 
BHES model. 

As for the physies eonelusions, the present paper updates and eorreets the results derived 
by the global BHES fit method whieh, following the eonsiderings just summarized, has been 
eompleted with an iteration proeedure in order to eaneel out possible biases. One thus eonfirms 
that almost all of the existing data samples eovering the annihilation ehannels with the 7 r° 7 , 77 , 
K^K~, final states and the dipion speetra in the —)■ 7 r=*= 7 r°z/ deeay aeeo- 

modate perfeetly the BHES framework. In the line of our previous works, one also finds that 
among the data samples eovering the —?• tt+tt” annihilation, the data samples provided 

by CMD2 and SND, the KEOEIO and now also the KEOE12 and BESSIII samples behave 
eonsistently with eaeh other and with the other eonsidered data eovering the various ehannels 
entering the BHES seope. 

The present update, whieh also ineludes the reeently published KEOE12 and BESSIII 
7 r+ 7 r“ samples, supersedes our previous results; these are mostly given in Table |3] and in Eqs. 
([TB . Erom a theoretieal point of view, it is interesting to note the eorreeted values for the c/s 
eoeffieients of the anomalous (EKTUY) terms of the HES model ifTSl [TSll : The eombinations 
c+ = (c 4 + C 3)/2 and ci — C 2 are found very elose to the usually assumed value, i.e. 1 ; in 
eontrast, c_ = (04 — c^)/2 = —0.166 ± 0.021 is non-zero with a Scr signifieanee. 

the case of a constant scale uncertainty, as for the CMD2, SND and BESSIII data, there is only one scale 
factor A. For most ISR data samples, the expression is slightly more complicated but easy to derive (see also the 
Appendix to ||35|) and the conclusions are obviously likewise. 
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Figure [3] displays the values for a^(7r7r, [0.63,0.958]) GeV derived from iterating the fits 
with the various available data samples. One observes a strong reduction of the uncertainty 
compared to the corresponding experimental value (about a factor of 2.5) and there is a close 
agreement between central values for all samples (or combinations of samples) which yield 
a good fit probability. The difference between the central values for the starting fit and the 
iterated one tends to indicate that biases are limited; this should be a consequence of also 
dealing with a large number of samples where the overall scale uncertainties are dominated by 
random statistical errors, as argued in the Appendix. 

Figure m exhibits the values for the muon Aa^ = when various combinations of 

7 r~^7r~ and samples are used in the iterated global fit procedure. The 

present study confirms that, within BHLS and because of its specific isospin breaking mech¬ 
anisms, one does not observe any serious mismatch between fits with only annihilation 
data and fits where these are supplemented with the r dipion spectra. The central valued for 
a^(e’''e“) and a^(e’''e“ + r) only differ by 2 units (NKS), 1 unit (NSK-i-KLOE-i-BESSIII-i-r) or 
0.7 unit in the global fit of all data samples (including BaBar) as can be seen in EigurelH 
Eigure |4] displays the value for Aa^ derived using all data samples except for KEOE08, 
which can be written : 

Aa^ = 37.02 + + [IHIftvsb ± 4.03*;, ± 6.3^.*,, 

where an estimate of the magnitude of possible uncertainties coming from outside the BHES 
framework is proposed. This exhibits a 5cr significance (which may reduce to 4.6(j - in the 
least favorable case - if the additional systematics are added linearly and assumed to play as a 
shift). One should note however that the fit probability is poor. 

The most probable value for the muon Aa^ is obtained by using the CMD2, SND, KEOEIO, 
KLOE12 and BESSIII samples - and the r spectra; this leads to : 

Aa^ = 37.68 + ± 4.12*;, ± 6 X.p. 

This BHES preferred estimate exhibits a 5.a significance for a non-zero Aa^, which may 
reduce to 4.7cr if one takes into account, as just above, the possible additional systematics. 
This solution is associated with a 99% fit probability. 

As a summary, even complemented with an iterative procedure shown in the Appendix to 
remove biases, the BHES approach favors a significance for Aa^ above the ~ 4.5cr level; this 
value is a lower bound obtained by including possible additional systematics added linearly. 
New data expected soon may further clarify the picture. The uncertainties now become sharply 
dominated by the region above 1.05 GeV, i.e. outside the BHES scope. 
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Appendix : Monte Carlo Tests of the Iterative Procedure 


A.l The Test Method 

In order to test the iterative method, one has developped a minimization eode whieh deals 
with speetra generated from a given underlying funetion Mtrue{s) where the parameters {a*} 
(whieh, of eourse, are known at the generation level) are fitted within the eode. The ’’experi¬ 
mental” speetra feeding this eode are generated using the true distribution smeared by introdue- 
ing gaussian uneertainty distributions. Indeed, for the purpose of testing our analysis method, 
it is eertainly the most appropriate to rely on ’’perfeet” data samples, with perfeetly known 
properties. 

For sake of simplieity, at the generation level, any ’’experimental” speetrum E is ehosen 
to earry 100 ’’measurements” mf, performed at 100 equally spaeed energy squared Sj points 
{si E [0, 1] GeV^), the same sequenee for all speetra. The ’’measurements” are derived by 
smearing the theoretieal values Mtrue{si) in the following way : For eaeh speetrum E, one 
assumes the ’’measurements” are sampled out from gaussian distributions in the following way : 

mf = Mtrue{si)[l ++ r]e'’fat{0A)] , i = 1, • • •, 100 (25) 

where 1) indieates the sampling on a gaussian distribution of 0 mean and unit stan¬ 

dard deviation generating the statistieal error; it varies independently from ’’measurement” to 
’’measurement” and from speetrum to speetrum. T]Mtrue{si) denotes the statistieal error eom- 
mon to all m*, g being some fixed fraetion of the order of a few pereents, ehosen the same for 
all the ’’measurements” in the speetrum E. 

On the other hand, Xe = 1) is the seale uneertainty affeeting speeifieally the 

speetrum E; as indieated by its definition, it is sampled out from a gaussian distribution of zero 
mean and a standard deviation. The overall seale uneertainty affeeting E is obtained via one 
sampling of 1) whieh, thus, earries the same value for all the ’’measurements” mf in 

the speetrum E. Of eourse, when going from a speetrum E to another E', another sampling 
of sfcaiei^i 1) should be performed. For speeifie tests, the overall seale uneertainty ean be 
switehed off {a = 0). 

One defines Nrep replieas (generally 1000) of Nexp (generally 5) experimental speetra eon- 
strueted as shown in Eq. (l25l) and submitted to a global fit where the parameters entering 
Mtruei^) j'JSt the parameters to be derived from the fit. The ’’true” statistieal error eovari- 
anee matrix Vij = [r]Mtrueisi)fSij is praetieally approximated by Vij = [gmffdij; we have 
avoided the unessential eomplieation of non-diagonal eovarianee matrix. The fit results de¬ 
rived for eaeh repliea are stored and then used to eonstruet the statistieal plots - true residuals 
and pulls -with the help of the known parameter ’’true” values. 

Therefore, we are just in the eonditions deseribed in Subseetion 14.2[ One should note 
that the MINUIT eode we have built performs the minimization of the samples and runs 
sequentially to treat the Nrep replieas within the same job. 

So, for eaeh repliea, the global minimized by our Monte Carlo MINUIT proeedure is 
simply a sum of Nexp terms like Eq. dH): 


X 


2 


— exp 

E xl 


(26) 
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When initializing the iteration proeedure, one uses Ae = rriE, i.e. the speetrum E serves to 
eonstruet its Xe; so differs from some other Ae> by statistieal fluetuations. When iterating, 
at first or higher order, they beeome identieal as Ae = Ae' = Mfu ~ M{afit). 

Obviously, eaeh sueh run provides simultaneously all the information allowing to examine 
the statistieal properties of the iterative method eorresponding to a given theoretieal ehoiee 
Mtrue{s)- The eomputer eode also allows an easy ehange of the funetional form of Mtrue{s) in 
order to examine the behavior of various kinds of non-linear parameter dependenees. 

The behavior of the fit parameters eompared to truth is, of eourse, the subjeet of the analysis; 
however, those of ’’physios quantities” derived from them are as important. For this purpose, 
we ehose to examine the ratioS : 

j- ^ Id Mfit{s)ds 

/o Mtrue{s)ds 

whioh has properties similar to those of the a^('Hi)’s, as the weighting faetor K{s) in Eq. ([I]) 
is an unessential eomplioation while looking for possible methodologioal biases of the iterative 
method. 

A.2 The Test Results 

The aim of the present Appendix is to report on numerieal analyses performed in various 
eonfigurations in order to examine how overall (global) normalization uneertainties and biases 
are related and whether non-linearities in the model parameters to be fitted lead to signifieant 
ineorreet estimates of errors. As Referenee PP71 whioh is faeed with the same kinds of issues 
as the present work, we do not plan to establish rigorously general theorems on these topios 
- assuming the soope of the issues would permit it. Nevertheless, one ean think that studying 
methods by relying on Monte Carlo teohnios is an aooeptable way to oheek its (praetieal) va¬ 
lidity under eommon oonditions. After all, the faot that Eq. ([3]) with A = M (the theoretieal 
funetion) is eonsidered free from biases is not weakened by the faet that the general (formal) 
proof of this property - if established - is not eommonly referred to. 

A.2.1 Analytical Shape of the True Distributions 

In order to use eonfidently fit results derived using the iterative method, one should examine 
the effeets of non-linear dependenees upon the fit parameters within eontexts similar to our 
physios distributions. The lineshape of the pion form faetor as a funetion of s on a given 
interval ean be qualitatively reproduoed using polynomials, ratios of polynomials, exponential 
of polynomials, sums of a Breit-Wigner funetion with polynomials ete ... with appropriate 
numerieal parameter values. 

We have applied the method outlined in Subseetion lA.il to perform fits relying on an inten¬ 
sive use of the tools provided by MINUIT taking various kinds of funetions Mtmeis), resem¬ 
bling - sometimes weakly - the pion form faetor. Running in sequenee migrad/Hesse and 
MINOS, we did not observe signifieant departures (beyond statistieal fluetuations) from equality 
between parabolie and MINOS errors; as the issue was to examine effeets of non-linear parame¬ 
ter dependenees this exereise was performed assuming statistieal uneertainties only. Therefore, 

^"^Remind that 0 and 1 GeV^ are the energy squared limits of the generated spectra. 
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this led us to conclude that, for the kind of experimental distributions one deals with, non-linear 
effects are not generally significant. For instance, using : 

Mtrue{s) = -xo u2 + c + ds + es^ , (28) 

[s — ay + ¥ 

T] = 3% and no scale uncertainty (to discard any need for iterating), the probability distribution 
was observed flat and the parameter pulls consistent with normal gaussians G{m = 0, cr = 1); 
the distribution of the ratio X for the 1000 replicas was also found well centered at 1 (actually 
its mean is 1.0001 and its standard deviation 1.62 x 10“^ from a gaussian fit with x^/^points = 
8.9/11). So, except for pathological cases which may always occur, non-linear dependences 
do not look practically an issue. 

From now on, we limit ourselves to reporting on using Mtrue{s) as given by Eq. (|2^ . 
Moreover, for sake of succinctness, we may only mention the fit parameter residual and pull 
distribution properties qualitatively and concentrate on discussing the distribution of the ratios 
X which, in fine carries - summarized - the relevant information. Each value of X entering 
this distribution is computed from a MINUIT fit of N^xp = 5 data samples and this is done for 
Nrep = 1000 replicas to construct numerically its distribution. 

A.2.2 Normalization Uncertainty and Iterative Method 

We first examined the results derived by fit of spectra with data points generated as in Eq. 
(1251) with a statistical uncertainty rj = 3% and generating the scale uncertainty A with a = 5%; 
so 7 ] is smaller than a. In this case, the interesting plots are gathered in Eigure[6l 

As one knows Mtrue{s), one can construct the N^xp partial with A = Mtrue(s) (see 
Eq.©) and minimize their sum using MINUIT. In this case, no bias is expected W2[ |45l |4^ 
and this is indeed confirmed by the top left panel in Eigure[^ where the distribution of the N^ep 
values for X is displayed. 

When, instead, one uses A = m (the data spectrum), the results are shown in the top right 
panel of Eigure[^ where one observes a shift of the central value by as large as 20% ! Denoting 
the result of the corresponding fit by Mq, one restarts fitting the same data by setting A = Mq, 
this - first - iteration leads to the distribution shown in the bottom left panel of Eigure[^ which 
looks identical to having used A = Mime- Denoting the fit solution of this first iteration by 
Ml, one restarts fitting the same data by setting A = Mi, and get the step 2 solution M 2 which 
correponds to the bottom right panel of Eigure which clearly indicates no change for the X 
distribution. 

So, one may conclude that the iterative procedure has already converged at the first iteration 
and so, we have Mi = Mtme- This fortunate high convergence speed has also been observed 
by [|471 and it is quite remarkable that this has allowed to recover frorr@ a 20% bias! 

Eit residuals are observed unbiased and pulls consistent with normal centered gaussians for 
A = Mtrue^ ^ = Mq and A = Mi. As for the yfi probability distributions, for A = m, it 
exhibits a huge spike at 1, while it is consistent with flatness (mean ~ 0.5 and r.m.s. ~ l/\/l2) 
for all the other cases. 

^^The numerical importance of this bias is intimately related with the ratio cr /77 = 5/3; if instead one works 
with cr /77 = 1, the bias coming out from fitting with A = m would only be 4%. 
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This already indicates that starting with A = m (the measured data spectrum) and iterating 
only once allows to give up knowing the theoretical function M beforehand to drop out bi¬ 
ases in physics quantity estimates. Moreover, as the parameter pulls are centered gaussians of 
unit standard deviations, the uncertainties derived from from the fit parameter error covariance 
matrix are reliable. 


FIT OF THF DISTRIBUTION INTFGRAL 






Figure 6: Distributions of the ratios X derived by varying the function A in the expression 
as indicated in each panel. The choice A = (i.e. the ’’measured” data sample) exhibits a 

20% bias while the other choices are unbiased. For more comments, see Subsection lA.2.2[ 


A.2.3 Effects of Subsamples Free from Normalization Uncertainties 

In the specific problem of globally fitting a large number of experimental data samples, one 
is faced with as many as 40 to 50 spectra to be treated [l34l [571 [5^ . Within this ensemble 

of data samples, one observes several configurations concerning uncertainties : some samples 
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have statistical errors dominated by scale uncertainties (the ISR collected data samples), while, 
in contrast, some others are reported with scale uncertainties marginal compared to statistical 
errors (the 7P data, for instance); sometimes, no specific information is reported 

concerning scale uncertainties, as for the the r dipion spectra ll49l[50l[5T]| . 


FIT OF THE DISTRIBUTION INTEORAL 






Figure 7: Effect of having 1 (top panels) or 2 (downmost panels) data sample(s) among the 
fitted Nexp = 5 samples simultaneously fitted. Left plots report on fitting with A = M (the 
truth), right plots on fitting with A = (the measured spectra); in the former case no bias is 
observed, in the latter case, the bias happens to be much limited. See text for more details. 


This makes interesting to examine configurations mixing samples of both kinds. In this 
paragraph, one summarizes the results obtained by running Nrep replicas of ensembles of 4 
data sets where, as before, the scale error is cr = 5% and the statistical error 7 = 3%, together 
with 1 data set with a = 0% (no scale uncertainty) and 7 = 6% (twice worse statistical 
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precision). This (4,1) combination will be supplemented with a (3,2) combination with the 
same characteristics. The main results are shown in Figure|71 Here we do not report on iterating 
the fit procedure, as obviously the results will follow the pattern shown in Figure [6l 





Figure 8: Probability distributions when fitting with A = Mtmth (left panels) or A = ruE (right 
panels). The top panel plots correspond to the case when among the N^xp = 5 fitted spectra, 
one is systematically free from normalization uncertainty; in the downmost panels 2 of the 5 
fitted spectra are free from normalization uncertainty. 


The top panels in Figure |7] display distributions of the ratios X in the (4,1) configuration. 
The left plot shows the case when the Nrep replicas are fitted using A = Mi^uth in the 
expressions. In this case, the absence of any bias is confirmed by the gaussian fit result shown 
within this plot. While using A = m®, the top right panel exhibits a 1.3% bias. Therefore, the 
effect of a single spectrum free from scale uncertainty out of 5 is enough to lessen dramatically 
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the observed bias : It reduces from 20% to 1.3%. 

The downmost panels in Figure |7] display the corresponding results when fitting N^ep repli¬ 
cas of (3,2) combinations. In this case, using A = in the minimized expression, leads 
to an even smaller bias (0.5%). 

So, even if they carry a poor statistical precision, having some spectra free from a (signif¬ 
icant) scale uncertainty is quite helpfull to strongly limit the real magnitude of a possible bias 
for a derived quantity. It is a quite interesting property to observe that some spectra with de¬ 
graded statistical quality supplementing other spectra dominated by scale uncertainties might 
be enough to avoid the need of an iteration procedure to unbias physics pieces of information. 

As for the probability distributions, comparing of the corresponding left and right panels in 
Figure [8] clearly shows that the departures from uniformity {i.e. average=0.5 and r.m.s.=0.289) 
due to using A = are quite limited. 

Nevertheless, when dealing with true experimental data (and thus unknown truth), one can¬ 
not take as granted that the number of samples with negligible scale uncertainties compared to 
statistical errors is sufficient to ascertain that biases are negligible. Therefore, in the practical 
case of the global fit of real experimental data performed within BHLS, secure results can only 
be ascertained by iterating until the change of a^{Hi) is small enough. 
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